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ABSTRACT 


The mathematical pendulum is one of the most widely studied problems in 
engineering physics. However, this is primarily limited to the classical pendulum 
with a single massless bar and mass configuration. Extensions to this include 
multi-degree of freedom systems, but many of the classical assumptions, such as 
a single bar per mass, are preserved. Several designs used in practice utilize 
multiple bars or trapezoidal configurations to enhance stability. Such designs 
have not been studied in detail. Also, there is a need for additional work to fully 
analyze their response characteristics. The two-string pendulum design 
characteristics are initially investigated, both in terms of oscillation and string 
tension. Analytical and numerical methodologies are applied to predict the 
response of the two-string pendulum in free and forced oscillations. Results are 
validated utilizing comparisons generated by simulations conducted with a 
standard commercial software package. A preliminary optimization study is 
conducted for a driven two-string pendulum. Finally, in a typical design case, it is 
shown how to apply the results of the analysis and optimization studies 
developed in this work. 
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I. 


INTRODUCTION 


A. BACKGROUND 

The pendulum is one of the most commonly studied mechanical 
structures. Galileo himself made the first observation of the pendulum’s 
performance. He set the roots of the pendulum’s physics by realizing that the 
frequency of oscillation for a given pendulum is independent of its displacement’s 
amplitude for small amplitude oscillations, which enables the linearization of the 
equations of motion. [1]. 

Following this initial observation, numerous studies were performed on the 
pendulum’s physics. These studies were not limited to the simple or standard 
pendulum, but were extended to more complex pendulum structures. Some 
examples are the bifilar pendulum, the torsion pendulum, the double pendulum, 
and the Foucault pendulum. 

To describe analytically the oscillatory behavior of the numerous types of 
pendula, several techniques were developed. These techniques attempted to 
predict the solution of the pendulum’s dynamics, either by treating its attitude as 
a linear approximation, which produced acceptable and fairly accurate results for 
a limited number of cases (linearized pendulum, , etc.), or by facing its nonlinear 
character recruiting and applying complex methods (perturbation methods or the 
Lagrange method, which produces the exact nonlinear equation of motion). 

From a mathematical analysis viewpoint, the more composite the structure 
of the pendulum, the more complex the solution approach method. Some types 
of pendula present a linear character when the oscillation is constrained in small 
displacement amplitudes (standard and double pendulum). Large amplitude 
oscillations are nonlinear and sometimes show chaotic behavior (double 
pendulum). If the problem is approached from an energy perspective, in linear or 
linearized undriven oscillations, energy is conserved. Nevertheless, there are 
examples where energy is not conserved. Such cases may be simple to analyze. 
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such as an oscillation with a constant damping coefficient. Others, such as 
oscillations, including jumps or variable accelerating support motion, are not 
simple to analyze. 

B. MOTIVATION - OBJECTIVE 

Several mechanical structures have dynamic behavior similar to that of a 
pendulum. Some examples are different types of cranes and crane-like 
mechanisms used for marine vehicles, such as AUVs, rescue boats, and naval 
vessels for launching and recovery. Their multi-rope suspension mechanism is 
not designed to oscillate; rather, it is designed to prevent swinging and to 
promote stability. Another example of a pendulum system designed to prevent 
swinging is the Blackburn pendulum. [1]. 


A 



Figure 1. The Blackburn pendulum 

The Blackburn pendulum’s suspension system is comprised of two strings 
forming a "Y” shape. The purpose of this design is to provide two distinct and 
effective pendulum string lengths. Therefore, for an oscillatory motion within a 
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plane defined by the Y shape, the length of the pendulum is the distance from 
point B to point C. In contrast, for a motion within a plane perpendicular to the Y 
shape, the effective length is the distance from point A to point C. Therefore, two 
distinct frequencies of oscillation define the motion of the Blackburn pendulum in 
the two perpendicular planes. Thus, its modes of oscillation are described by 
superposition of the two basic frequencies. 

Nevertheless, all the previously mentioned characteristics of the Blackburn 
pendulum depend on the assumption that point B will not be able to oscillate with 
respect to one of the two strings that starts from B and ends at the support 
points, i.e., the distance of B from the two support points will always be constant. 
However, this assumption is valid under restrictions. A violent motion of the 
support would theoretically be able to disturb the partial equilibrium of the "Y” 
pendulum. 

The motivation for this thesis was to obtain analytical expressions, as well 
as numerical solutions, to describe the behavior of oscillating pendula with 
designs similar to the one of the "Y” pendulum. The basic model design that will 
be initially analyzed is the one of a Blackburn pendulum in which point A 
coincides with B, and is allowed to oscillate only within the plane defined by this 
“V” shape (more detailed description will be provided in the following sections). 
Moreover, since such a "V” design is the basic element component of more 
complex suspension designs, its capability of preventing or damping oscillations 
will be investigated. 
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THE TWO-STRING PENDULUM 


A. MATHEMATICAL DESCRIPTION OF THE TWO-STRING PENDULUM 

1. Introduction 

The purpose of this section is to analyze the motion equations of the two- 
string pendulum. For brevity within this thesis, this will be referred to as the 2- 
pendulum. Several techniques will be used to analyze the motion of the 2- 
pendulum with the aim to present an accurate solution to the nonlinear equations 
of motion that describes its motion. The main purpose is to compare main 
features of the dynamic response of the 2-pendulum with those of the standard 
single-mass, single-string pendulum. These comparisons will be in terms of the 
amplitude of the resulting motion, the resulting natural frequency of the oscillator, 
as well as the string tension. Different values of initial velocity and initial angular 
displacement will be used to investigate the system’s response to both an initial 
excitation and an initial displacement from its equilibrium position. 


2. The Standard Pendulum 



Figure 2. Small angle oscillation of the standard pendulum 
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First, we will investigate the motion of the simple gravity pendulum 
undergoing free undamped vibration. The model that we will use is shown in 
Figure 1 -- a sphere of mass m which is suspended with a massless string of 
length = t from its center of gravity. 

To derive the equations of motion, we will use the moment method. 
According to Newton’s Law, the inertial moment of the mass m must be equal to 
the sum of moments applied on the mass. 

=-mg£ sin 0 (2.1) 

mP0 + mglsm0 = O (2.2) 

Equation (2) above is nonlinear since it contains the term sin^. We can 
substitute the trigonometric nonlinearity with its Taylor series expansion 

n5 

sm0 = 0-^+^... (2.3) 

To simplify our calculations, we will initially assume only small amplitude 
oscillations. In this case, we can assume without compromising accuracy that 

sin 6'-6' (2.4) 
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Figure 3. 


Linear-nonlinear discrepancy 


As demonstrated by Figure 3, there is a range of 6 where equation (4) is 
valid and provides a good approximation. As depicted in the same figure, the 
difference between the linearized and the nonlinear expression is negligible for 


values of 


71 
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For such small amplitude oscillations, the equation of motion can be 
linearized 


or 


+ mgi6 = 0 


0+^e=Q 


(2.5) 

( 2 . 6 ) 


The solution to the above differential equation for initial displacement and 
velocity of 
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becomes 


6»(0)= e, 

m =4 


(2.7) 


e= 


m 

(On 


sin coj + ^(0) cos coj 


( 2 . 8 ) 


where 



(2.9) 


is the natural frequency of the linear pendulum. 
From Figure 2, the restoring moment is 

= mg sin 6 


or, for small oscillations, 


= mg9 


( 2 . 10 ) 


( 2 . 11 ) 


Another feature of the pendulum worth mentioning is the pendulum jump 
[5]. Pendulum jumps occur when the mass leaves the circular path defined by the 
string and follows other trajectories. Such response is only observed if the string 
tension becomes zero. This is dependent on the position of the pendulum and 
the initial conditions. From Figure 2, we can visualize that the centripetal force, 
which obliges the mass to follow a circular path, is 

2 2 

j-i TTtli rj-f ^ ^ ^ rj-f TTlli ^ ^ 

Fc = ^— = T-mgco^U <^1 =—^ + m^cos6' (2.12) 


where u is the linear velocity of the mass. The magnitude of u depends only on 
the initial conditions, since only conservative forces govern the pendulum’s 
motion. Assuming the initial conditions (2.7), which are equivalent to 

Zq =.^(1-cos^q) and Uq=Wq, the initial energy of the mass is 


EQ=TQ+VQ=]^mul- 


■mgZQ- During oscillation, the energy of the mass for an 
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I /y 

arbitrary angle 6 would be E = T+V =—mu +mgz , where z = i(l-cos0). 

Since the total energy is conserved throughout the oscillation (absence of 
damping or other non-conservative forces) £" = , the magnitude of the linear 

velocity would be 

u = 2.^-gz (2.13) 

V m 


Finally, the tension of the string throughout the oscillation, from (2.12) and 
(2.13), will be equal to 

2(EQ-mgz) 


T = 


£ 


+ mg cos 0 


(2.14) 


where a positive sign denotes orientation towards the center of oscillation. This 
promotes the circular path motion of the mass. The only term of the latest 
equation that may take negative values and, thus, lead to zeroing the string 
tension, ism^cos^. It may only happen for angular displacements greater than 


n 

2 ■ 

The possibility of jumps is present only if the initial conditions allow the 

oscillator to climb to After that point, the path of the mass would be 

governed by the regular equations which define a projectile with a given initial 
linear velocity magnitude and orientation; that is, until the position of the mass 
will again reach a distance equal to the natural length of the pendulum’s string. 
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3. 


The 2-Pendulum 



Figure 4. The 2-pendulum 


The geometry of the 2-pendulum is pictured in Figure 4. The main 
distinguishable characteristic of the 2-pendulum, compared to the simple 
pendulum, is that the mass m is now suspended by two strings. To be consistent 
with the previous model, the strings’ perpendicular distance from the plane of 

£ 

suspension will be /. Consequently, 1 = -. The point of application of the 

cos a 

strings on the mass is again at its center of gravity while each string forms an 
angle a with respect to the perpendicular. Our aim is to evaluate the effect of this 
angle on the motion of the mass m. 

For the purposes of this preliminary analysis, we will assume that the 
governing force during the pendulum oscillation is the weight. That means that 
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we will not, for now examine the forces applied on the mass to the point where it 
reaches the lower point of its trajectory and switches the string about which it 
swings. 



gcos{6+a 


Figure 5. Small angle oscillation of the 2-pendulum, @<0 
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Figure 6. Small angle oscillation of the 2-pendulum, 0>O 


Assuming an angle 6 (counterclockwise for positive values), as in Figure 
6, and equating inertia and non-inertia forces applied on the mass as before, the 
resulting equation of motion is 

ml^6 + mglsm(^0 + a^ = O (2.15) 

or 

d + Y^in(0 + a) = O (2.16) 

To further develop equation (2.16), we should take into account that the 
term ^m{6+a) takes positive values for the case of 6>0 (Figure 6), and 

negative values for the case of 6<0 (Figure 5). Flowever, since the restoring force 
should always oppose further increase in the magnitude of the angle 0, it takes 
the same magnitude for the same||a-H^||. Furthermore, using the trigonometric 

identity sin(a-hZ?) = sinacosZ?-hcosasinZ?, equation (2.16) becomes 


^-hy(sin6'cosa-hcos6'sina) = 0 


(2.17) 
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Examining equation (2.17), we can see that the consistency of the 
restoring force orientation is preserved in the term sin^coscr (because of 
sin^), but not in sinacos^ (since cos0 does not change sign when ©changes 
sign). Consequently, to preserve the restoring force consistency for all values of 
6, we should rewrite equation (2.17) as follows 


^ + 4^sin6'cosa±4^cos6'sina = 0 


This is equivalent to 


6^ + y sin 6cos a + y cos ©'sin a sign 6 = 0 
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(2.18) 


The sign nonlinearity assumes the value (+1) for positive and (-1) for 
negative arguments. In addition to the above, if we assume small angles of 6, we 
can rewrite (2.18) as follows 

^ + y6'cosa + ysinasign6' = 0 (2.19) 

The restoring force from equation (2.15) is 

= mg sm(6 + a) (2.20) 

For small oscillations, equation (2.20) becomes 

= mg6cos a + mg sin a sign 6 (2.21) 

Our goal is to solve analytically and numerically differential equations 

(2.18). 

Finally, based on the previous analysis performed on the jumps for the 
standard pendulum, we realize that the 2-pendulum may exhibit the same 

TT TC 

behavior, only now the required angular displacement \%9+a> — 9 > —-a . 


4. Comments 

Despite its similarity with the standard pendulum, equation (2.18) exhibits 
some major characteristic differences. The most significant difference between 
the standard pendulum and the 2-pendulum is that the latter possesses a 
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discontinuity in its stiffness around the equilibrium position. As can be seen from 
(2.19), there is an instantaneous jump in the value of the restoring moment 
around zero. This is depicted in Figure 7 where the normalized spring constant of 
the 2-pendulum is shown for different values of a , along with the linear spring 
constant of the standard pendulum. This discontinuity at zero prevents us from a 
classical linearization of the 2-pendulum and requires different techniques for the 
analysis. Such techniques, based on the perturbation theory, are the subject of 
the analysis that follows. 

1.5 

1 

0.5 

0 

- 0.5 

-1 

- 1.5 

-1 - 0.8 - 0.6 - 0.4 - 0.2 0 0.2 0.4 0.6 0.8 1 

Figure 7. Normalized spring constant of the 2-pendulum 
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5. 


Velocity and Force Vector Analysis 



Figure 8. General force and velocity diagram at the switch point 

Before advancing with the solution of equations (2.18) and (2.19), we 
should take into account the dynamics of the mass m at the point when it 
switches strings about which it swings. In Figure 8, we can visualize the forces, 
weight B, instantaneous string tension Tinst (blue with continuous line vectors), 
and their resultant R (blue dashed line vector). During the change, it acts on the 
mass and actually contributes to not only the switch at the center of rotation, but 

to the velocities just before Vb (red vector), and just after, Va (green vector) -- 
the switch. Nevertheless, since the weight of the mass is a conservative force 
and does not undergo changes in magnitude during this switch, we can safely 
assume that its contribution in the velocity direction switch is negligible. 
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Figure 9. 


Force and velocity diagram at the switch point for a=45° 


According to the vector diagram in Figure 9, there is a limiting value of 
angle a, which, if exceeded, the string switch is impossible. In the instance that 

Tjnst lies in the same direction, but in opposite orientation to the Vb, this force is 


unable to change the direction of the velocity from vb to Va , irrespective of its 
magnitude. This limiting value is shown to be a=45°. In this case, this exact value 
of the mass will either stop or bounce in the same direction. However, the 


orientation of Tinst, will have a velocity magnitude less than or equal to Vb 
(depending on the energy absorption of the string). As a result, the following 
preliminary analysis. Chapter II, Section A, numbers 6, 7, and 8, will be valid only 
for a<45°. Furthermore, we will not consider the exact forces that affect dynamics 
and energy equivalence at the switch point. Rather, we will only take into account 


that their effect is the jump of the velocity from Vb to Va , where 


Va 



6. Numerical Integration Solution 

The numerical solution is performed using the Matlab function ode4 5. 
This ode4 5 Matlab function solves numerically the equation of motion and is able 
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to take into account nonlinearities, such as trigonometric functions and sign 
changes. Thus, ode4 5 will provide the exact solution, provided that an 
appropriate, fine enough, time step is selected. The Program_1.m Matlab-file 
was used to plot the angular displacement 6 versus time for different values of 

angle a, assuming initial conditions 0{Qi) = 0q and ^(0) = 0. In parallel, we 

plotted the angular displacement 6 versus time for a simple pendulum of string 
length £. In addition, the Program_2.m provided the same kind of plots, where 
the initial conditions were ^(0) = 0 and ^(0) = 0q . 


7. Analytical Solution 

Three different methods will be used to attempt to derive the exact or 
approximate solution to equation (14) or its equivalent (16). 


a. Exact Solution for a Time Increment 

One method we could use to obtain the solution of equation (2.18) 

is to solve it for the time period where the term y sin crsign ^ has a constant 

sign [2]. Thus, assuming an initial angular displacement 6(0)=A, the equation of 
motion for O<0<A is 


or 


^sma 

0 = --^ -+ 




8_ 

I 


cos a 


A + 


ySincr 

^coscr 
/ y 


cos 


lycoscr 




( 2 . 22 ) 


^ = - tan cr + (A + tan cr) cos 


-^cosa 


I 


(2.23) 


Let us name the square of the pseudo-natural frequency of 
oscillation and of the equivalent pseudo-forcing amplitude as bellow 


cDq =yCosa 
g . 

P = ySincr 


(2.24) 
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At the time t = — arccos 


CO, 


0 


r \ 

1 

11 M 

V P J 


when 6=0, the velocity is 




^4 

A +—Y 
V 


sin co^t 


and then the equation of motion becomes 


with initial conditions u{t^^ = -coQ 


e+cole=p 
\ P ^ 

A H- IT 

V 


(2.25) 

sin<z>of and 6{t^) = Q, and so on. The 


value is equivalent to one quarter of the period of oscillation. 
Finally the period of oscillation is 


T = —arccos 
exact Q) 


0 


r \ 

1 

1+^ 

P j 


or 


exact 


ycoscr 


r arccos 


r ^ 

1 

V tan a y* 


so the frequency of oscillation is 


CO 


Itt 


exact 7 " 


exact 


(2.26) 


(2.27) 


(2.28) 


b. Method of Slowly Changing Phase and Amplitude 

A second method for obtaining the analytic solution only for the 
linearized equation (2.19) or, considering (2.24) 

e + co^0+Psign0 = O, (2.29) 
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is using the method of Slowly Changing Phase and Amplitude [2], This method 
states that the equations of the slowly changing phase y/ and amplitude a are 


i = -—^^-Psign |a sin (6'+^)} J sin (a + ^) tia 


(2.30) 

(2.31) 


Evaluating the above integrals, 


a = 0 and xj/= 


a = A and y/= 


AP 


27T2LCOn 


IP 


JlACOr 


■ + (^0 


n 


Finally, if we choose the initial conditions, such as the approximate 

solution is 


6^ = Acos 

and the approximate period is 


IP 

(^Q+ -T- 

nAco, 


oj 


Itt 


" app 


co^ 




1 + 


IP 

nco^A 
0 J 


(2.32) 


(2.33) 


while the approximate angular frequency is 


^app ^0 


1 + 


IP 

nco^Aj 


considering (2.24) 


G)app=.\^O^CL 


( 


2tancr 




(2.34) 


(2.35) 


c. Energy Method 

The application of this method suggests that we could approximate 
the equation of motion of the enhanced pendulum with one of a simple 

pendulum. However, the restoring force would be increased by a factor so that 
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the total work produced by this force along a time period would be the same as 
the forces that govern the motion of the enhanced pendulum. 

From equation (2.20), the work produced by the restoring force of 
the enhanced pendulum fora motion 0<6'<6'o is equal to 

\^^mg^m[a + 0)d0 (2.36) 

while the one of the simple pendulum is 

\l^mgm\0d9 (2.37) 

Dividing (2.36) by (2.37), we obtain a factor of 

F = mg + 0)d0!mg ^mOdO (2.38) 

which provides us with an order of increased magnitude offerees applied during 
the oscillation of the enhanced pendulum compared to the ones applied on the 
simple pendulum. 

Thus, we will multiply the term mg^mO by F and solve the 
approximate equation of motion 

|9 + F^sin6' = 0 (2.39) 

or in linear form, 

0 + F^e = {) (2.40) 

Thus, the approximate angular frequency of oscillation is 


CO, 


app 



(2.41) 


8. Evaluation of Results 

To evaluate both numerical and analytical approximations, we will first plot 
the three angular displacement solutions (ode45 numerical integration - Slowly 
Changing Phase and Amplitude method - Energy method). This will be 
compared to the solution of the standard pendulum in terms of gain and 
frequency of oscillation. This will best show that the 2-pendulum offers significant 
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gains in terms of increased restoring moment and frequency of oscillation. It, 
also, measures how close it comes to the exact solution (represented by the 
ode45 solution). The predicted response will be presented for two distinct initial 
conditions -- initial angular displacement and initial angular velocity. In addition to 
the above, we will plot the phase trajectories for the same initial conditions, as 
well as the comparison of the predicted angular frequencies of oscillation. This 
will help to better visualize the discrepancy of the prediction methods. 
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Figure 10. a=TT/6 
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Figure 11. a=n/12 



Figure 12. a=TT/18 
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Figure 13. a=TT/36 
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Figure 14. a=TT/90 
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Figure 15. a=n/180 
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Figure 16. Phase diagram a=TT/6 
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Figure 17. 


Phase diagram a=TT/12 
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Figure 18. Phase diagram a=TT/18 
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Figure 19. Phase diagram a=TT/36 
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Figure 20. Phase diagram a=TT/90 
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Figure 21 . Phase diagram a=TT/1 80 
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Figure 22. a=n/6 with initial velocity 
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Figure 23. a=TT/12 with initial velocity 
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Figure 24. a=n/18 with initial velocity 
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Figure 25. a=n/36 with initial velocity 
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Figure 27. a=n/180 with initial velocity 
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Figure 28. 


0o=Tr/6 
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Figure 29. 0o=Tr/1O 
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Figure 30. 0o=Tr/2O 
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Figure 31. 0o=Tr/3O 


In Figures 28, 29, 30, and 31, we can evaluate the accuracy of the 
response prediction method with respect to the initial displacement do and value 
of the characteristic angle a. According to equations (2.28), (2.35), and (2.41), 
the oscillation angular frequency depends on both parameters. Therefore, the 
four plots consider four arbitrary values for the initial angular displacements. This 
allows us to better appreciate the accuracy of the two prediction methods, SCPA 
and Energy, with respect to the main attribute of the 2-pendulum, which is angle 
a. We should note that the three equations are valid only for the linear 
approximations of the equations of motion. In addition, the graphs consider 
a^45°, a restriction we already mentioned in the beginning of the analysis. 

According to the plots, we conclude the following: 


32 























• The angular frequency w of the 2-pendulum is greater compared to 
the frequency of the standard pendulum which has the same equivalent string 
length i. . 

• A smaller maximum angular displacement for the 2-pendulum 
compared to the standard pendulum results from applying the same initial 
angular velocity. 

• The discrepancy between the two predictions, and the exact 
calculations, is very low for small values of a, but increases for larger ones. More 
specifically, while the Energy method follows the exact solution slope trend 
(although providing higher predicted values), the SCPA method increases almost 
linearly. 

• The SCPA method tends to provide more accurate predictions for 
larger values of do. If do becomes too small, the discrepancy from the exact 
solution grows and, therefore, the predictions are incorrect. 

9. Horizontal and Vertical Displacements 

Despite the obvious advantage that the 2-pendulum possesses in terms of 
increased angular frequency of oscillation, increased restoring force, and 
reduced maximum angular displacements, it is worth mentioning a disadvantage; 
the differences of both horizontal and vertical displacements associated with 
angular displacements. The following figure depicts the relationship that governs 
the maximum horizontal and vertical displacements for 6o=Tr/18 and a range of a. 
This figure depicts that, for each angular displacement of the 2-pendulum, the 
horizontal and vertical displacements are smaller and greater respectively when 
compared to the standard pendulum. This is due to the two pendula’s different 
reference systems used to measure their angular displacements. Therefore, to 
obtain a better appreciation of the comparison between the standard and the 2- 
pendulum, the previous examples of a will be plotted in terms of linear 
displacements instead of angular displacements (Figure 33). 
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Figure 32. Variation of maximum vertical-horizontal displacements wrt 

characteristic angle a 
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Figure 33. 2-pendulum’s mass trajectories during oscillation 
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B. STRING TENSION INVESTIGATION 


1. Analytical Formulation 

Another important concept that requires further investigation is the tension 
applied on the string of the 2-pendulum. As already seen in equations (2.15) and 
(2.2), the restoring moment on the mass m is larger compared to the standard 
pendulum. Contrarily, the string tension has smaller values since two strings, 
instead of one, are used to suspend the mass. As pictured in Figure 3, while in 
equilibrium, positioning the string tension is dependant only on the value of the 
characteristic angle a 


Ts = ^^ 
2cosa 


(2.42) 


According to Figure 4, during swinging and while the sign of angle 6 is 
continuous, the string tension can be computed as follows 


Ts = Fc + jmgcos{a + 0)j (2.43) 

2 

where Fc = —^— = m6^l is the centripetal force. 

Equation (2.43) shows that the magnitude of static Ts ranges from 

IflQ Jl 

<r5'<oo,for 0<a + 6< — , proving that increasing the characteristic angle 

does not have only positive results in the overall design. 

To determine the dynamic string tension, trigonometric equalities are 
applied and 

Ts = mOH + \mg cos a cos 0-mg sin a sin 6\ 

71 

To attempt to linearize the above equation for d< —, 

10 

Ts = m6H + \mg co^asignO-mg^ma6\ (2.44) 
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For the standard pendulum, the calculations are easier. Therefore, the 
tension applied on the string, while in equilibrium position, is 


II 

(2.45) 

and during the oscillation is 


Ts = mO^l + \\mg cos 6\\ 

(2.46) 

or, for small amplitude oscillations. 


Ts = md^l + mg 

(2.47) 


Since ode45 solves the exact differential equation of motion and, we 
already evaluated the accuracy of the predictions provided by the SCPA and 
Energy methods, the string tension simulation will be performed using the results 
provided by the ode45 method. 



Figure 34. Tension for a=TT/18 - 6o=Tr/10 
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Figure 35. Maximum and minimum tensions 



Figure 36. Maximum and minimum tensions 3-D variation 
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Figure 37. Minimum tension variation 



Figure 38. Maximum tension variation 
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In Figure 34, we should note that this is not the variation of tension in one 
of the strings, but, rather, the sum of them. Although, each string’s tension is zero 
for every other half period, the fact that we describe the motion of the mass as an 
equivalent standard pendulum, with altered restoring force, and we do not 
distinguish the angular displacements with respect to the different strings from 
each other, leads us to treat the string tension with the same methodology. As a 
result, in each peak in Figure 34, the mass switches strings, and these exact 
switches take place when the string tension takes its maximum value. 

From Figures 35-38, we can evaluate the variation of maximum and 
minimum tensions exerted on the strings with respect to both a and Go. On the 
one hand, as expected, increasing do increases the Tmax and decreases Tmin. 
On the other hand, increasing a, and keeping do constant, steadily decreases 
Tmin. An interesting point is that Tmax does not exhibit the same behavior. In 
Figure 39, we can see for which value of a we will have the maximum tension for 
a range of initial displacements. 
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Figure 39. Characteristic angle producing the maximum tension for given do 
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C. NUMERICAL SIMULATION 

A very precise method to simulate the response of the 2-pendulum is with 
the use of the VISUAL NASTRAN 4-D software. In this section we attempt to 
validate the methods that we used to predict the 2-pendulum response in 
Chapter II, Section A, by comparing their results with those produced by VN 4-D. 

1. VISUAL NASTRAN Software - SOLIDWORKS 

To simulate the response of the 2-string pendulum in real time, a specific 
procedure had to be followed. This procedure was initiated by designing a model 
of a 2-string pendulum using SOLIDWORKS. This 3-D design software allows 
the exact geometrical design of structural models, which are compatible with the 
VN 4-D software. This model is then transferred into VN 4-D. The distinct 
attribute of VN 4-D is to solve differential equations in 3-D space and time. The 
final model used in VN 4-D was comprised of geometrical parts: mass properties 
- developed in SOLIDWORKS -, and their assembly - developed partially in 
SOLIDWORKS and partially in VN 4-D. The assembly constraints, which resulted 
in the rigidity of the platform where the 2-pendulum support points were 
positioned, were developed in SOLIDWORKS. The string constraints were 
developed in VN 4-D. Before initiating each VN 4-D run, the mass had to be 
positioned in the exact coordinates. This was the equivalent to a combination of 

initial angular displacement 0^ and initial angular velocity 9^ ofthe VN 4-D. The 

results of the simulations (in terms of mass position wrt which is the model’s 
global system of axis, linear velocity, and tension of the strings) were then 
exported into Microsoft Excel spreadsheets. In addition, as defined in Figures 5 
and 6, the mass position results were processed to obtain the angular 
displacement of the pendulum. Finally, all the above results were graphed with 
the aid of Matlab. This allowed comparing them with our analysis results. 
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Figure 40. Front view of 2-pendulum VN 4-D model 
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Figure 41. Side and bottom view of 2-pendulum VN 4-D model 

The front, side, and bottom views of the model that were used for the VN 
4-D simulations can be visualized in Figures 40 and 41. In this model, one can 
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adjust the mass of the mass (blue sphere) and the characteristic angle by 
altering the length of the strings (black lines) or the distance between the support 
points (red cubes). In addition, an initial displacement can be applied to the mass 
by precisely positioning the mass to a specific position in terms of an exact set of 
coordinates (X,Y and Z), with respect to the Cartesian global system of axis 
represented by the black arrows. It should be noted that the exact positioning of 
the mass is not possible. For a few of the first time steps of the simulation, this 
results in a slight deviation of the expected results. The numerical integration 
parameters that were used in all the VN 4-D simulations are: 


A/A 

Simulation Parameter 

Value 

1 

Integration Method 

Kutta-Menson 

2 

Animation Time Frame 

0.001 sec 

3 

Integration Step 

Variable - 0.0005 sec 

4 

Steps per Frame 

2 


Table 1. VN 4-D simulation parameters 


The above parameter combination yielded results that were tabulated in a 
time-step size of 0.001 sec. This way compatibility was achieved with the 
previously presented results using the developed Matlab code. 

In the following example, we chose to solve the case of the 2-string 

71 

pendulum of characteristic angle «=— and initial angular displacement 


6'n 


n 

10 


. For consistency we used t = while the actual string length was 


/ = 1.0154m. 


The VN 4-D results will be compared to the predictions of the three 
methods described in Chapter II, Section A, numbers 6 and 7, for the same 
values I, I ,a and do. 
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Figure 42. VN 4-D - a=TT/18 - eo=TT/10 
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Figure 43. Theoretical predictions for a=n/18 and 6o=Tr/10 
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Figure 44. Angular frequency comparison for a =tt/1 8 and 6o=tt/1 0 
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Figure 45. Typical velocity variation for 2-pendulum 


2. Evaluation of Results 

In Figures 42 and 43, we can see the oscillation prediction using VN 4-D 
and the analytical methods, respectively. To more precisely appreciate the 
discrepancies between those methods’ predictions. Figure 44 will be used. We 
can clearly distinguish from Figure 44 that VN 4-D predicts that, with time, the 
amplitude of oscillation decreases and the angular frequency increases. Further, 
for angular frequency, the comparison was indispensable. We decided to plot the 
angular frequency every half-cycle to be able to capture the decrease in w in the 
most possible detail. From Figure 44, we can realize that ode45 provided the 
most accurate prediction for the w for the first half-cycle, compared to SCPA and 
Energy methods. Nevertheless, after that point, the w of VN 4-D started to 
increase almost linearly. This increase in w was not captured by the three other 
methods. 
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To understand why this increase in w takes place, we plotted the angular 
frequency of oscillation predicted by VN 4-D. As seen in Figure 45, for every half¬ 
cycle when the mass reached the point where it switched strings about which it 
oscillated, the velocity magnitude dropped by a certain amount. This 
discontinuity, just before and after the switch point in velocity magnitude, raises 
questions about energy conservation. In the following sections, we are going to 
investigate in more details the dynamics that surface at this point and their effect 
in terms of energy loss of the oscillation mass. 

D. THE REAL 2-PENDULUM 

1. Velocity and Force Vector Analysis-Impact Dynamics 



Figure 46. Velocity and force vectors for 2-pendulum at switch point 

In Figure 46, we can visualize the forces that act on the mass during this 
change and contribute to the switch of center of rotation, as already depicted in 
Chapter II, Section A, number 1 and Figure 7. Nevertheless, in this section we 
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are going to emphasize the way that the forces change the direction of the linear 

velocity from Vb \.oVa, as well as other effects, such as energy loss and 
termination of oscillation for a>45°. 

According to Figure 46, the dynamics of our problem strongly resemble 
the dynamics of a ball colliding with a rigid wall. For a<45°, the only known 

parameters are: i) vector Vb (direction, orientation, and magnitude), ii) the 

direction and orientation of vector Va, iii) the weight of the mass, and, finally, iv) 
the direction and orientation of Tmst- Based on these parameters, we will attempt 
to predict the response of the real 2-string pendulum. 

2. The Real 2-Pendulum with a<45 

In our initial analysis, we are going to assume that the string is stiff enough 
and allows no elongations even under large values of tension forces. Therefore, 
we are not considering displacements along the direction of the string. 
Nevertheless, we are going to assume that an amount of energy is lost during 
this switch of velocity direction. Following the discussion of Chapter II, Section A, 
for this case the pendulum will oscillate ideally in a manner that was already 
analyzed in that section. In addition to the previous analysis, we are going to take 
into account the energy losses that occur at the switch point due to the 
introduction of nonlinearity. For this purpose, we will apply the equation of motion 
of oscillation between two arbitrary consecutive switch point impacts, i and i + l, 

which is equation (2.18), ^-hysin6'cosa-hycos6'sinasign6' = 0, and the 

simple impact rule at the switch point i , for instance. The simple impact rule 
states that the impacts are instantaneous and the energy loss would be 
represented by a reduction in the velocity magnitude. That is, 

= l■0_, (2.48) 

where r is the coefficient of restitution at impact [3]. 
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Referring to Figure 46, we could be more specific and precise while 
calculating the energy absorption during impacts described by equation (2.48). 
As already mentioned, the impact exerted on the mass by the string can be 
described by a collision of a rigid ball with a rigid wall. Since the weight of the 
mass is a conservative force, it does not absorb any energy during the collision. 
Therefore, T/nsf is the only force that contributes to energy absorption, which 
equals the work done by force W during the period of contact [4]. Moreover, the 
sole velocity components that are affected by W are those parallel to the 
direction of Tinst. Dealing this time with linear velocities (instead of angular), we 

can see that the components of Vb and Va , which are parallel to the direction of 
the string causing the impact, (heavy dotted line in Figure 46), are , with a 


magnitude of 


Vb 


ajf 


Vb 


cos(90-2a), and 


Va 


Cijf 


Va 


cos(90-2«). Finally, using 


yiy = y-i.ya^^=V^j^+i and Vb^j^=V^j^+i 


Va^jj, with a magnitude of 
the notation Va=V_^-, 


where 




2 


(2.49) 


^inst=Te-i-T,,i =^ = Pe-i-Pc-i ( 2 - 50 ) 

and At is the infinitesimally small time duration of the impact. 

In equation (2.50), is the impact force exerted by the mass m on the 

string causing instant elongation. is the impact force exerted by the string 

contracting and releasing its strain energy on the mass m. are the 

respective impulses. Consequently, in the circumstance where there 

would be no energy loss during impact and we would deal with a completely 
elastic collision. 
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In addition to the above, from the point of view of the conservation of 
momentum and considering, once more, only the velocity components that are 
parallel to the string direction, during the phase of the instantaneous elongation 
of the string 


MV,=MV^^_. + P^, (2.51) 

where i^o is the velocity of the mass at the end of the elongation phase, is 

the reaction impulse force that is exerted on the mass, and M is the effective 

TYITY^ 

mass which equals M =- r. In this case, m! is the fictitious mass of the 

m + m 

object onto which the mass m collides. Since this object is represented by a rigid 
wall whose momentum is infinite, we assume that m=co. Therefore, M =m 
and, since Vq = 0, the magnitude of the impulse force is 

Pel=^^aff-i 

or 

p^i =mV_iCOs{90-2a) = T^_^At (2.52) 

where At is the infinitesimally small time duration of the impact [3]. 

The next step should be to calculate the pc = Tc + lAt. If we conduct a 
similar analysis, we would see result 

p^ =mV+,cos(90-2a) = r^^,Af (2.53) 

Looking more closely at equation (2.53), we realize that if after 

the string-change point, there would be a linear velocity component parallel to the 
string’s direction in addition to the linear velocity component perpendicular to the 
new string’s direction. The second component promotes the smooth circular 
motion of the mass, while the first component would push the mass to deviate 
from the circular orbit. Consequently, if we want to predict the path of the mass 
that would not deviate from the regular pendulum’s oscillation path we will 
assume that 
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(2.54) 


V^,cos(90-2a) = 0 

Thus, equation (38) should be manipulated and corrected as follows 

A _ _ j(V^^ sin(90-2a))^ + (V^- cos(90-2a))^ _ 


(V_^ sin(90 - 2a)y + (ry_, cos(90 - 2a))^ _ 

^^(sin(90-2a))^ + (rcos(90-2a))^ 

i9+, = V(sin(90- 2a))2 + (r cos(90 - 2a)f (2.55) 

where, considering (2.55), the coefficient of restitution should be r = 0. 

The cancellation of the parallel component of the velocity during impact is 
not only an assumption that facilitates the work of predicting the response of a 
real 2-pendulum, but, also, a very good approximation. The sudden impulse 
shock that is applied on the string from the mass produces a compressive wave 
that propagates throughout the length of the string. This wave will finally die out 
under the effect of damping. Damping in almost inextensible strings may be 
caused by several factors, such as viscosity of the air and internal friction. The 
velocity of elastic wave propagation within the string mass is similar to the speed 
of sound [5][6]. Therefore, by comparing this velocity with the time scale 
associated with the periodic motion of the pendulum, the time required for the 
completion of the energy loss related to this shock is negligible. All the above 
support the validity of the assumption that the coefficient of restitution, associated 
with the impact of the mass on the string, is zero. 

Finally, from equations (2.54) and (2.55), the amount of energy loss at the 
switch point is 

^loss =^mv^^_^((l-r)cos(90-2a)f (2.56) 
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The system of equations (2.18) and (2.48) or (2.55) can be solved 
numerically using a constant time step of A? = 0.001sec and the forward 
difference scheme. This takes into account equation (2.48) or (2.55). The Matlab 
code, Program_3.m, provides the solution to the system of equations for 
predetermined number of periods. 

In the following sections, we are going to compare the prediction of the 2- 
pendulum response in terms of angular displacement, angular velocity, and string 
tension. We will use the above m-file (solving equations (2.18) and (2.55)) 
combined with those obtained from VISUAL NASTRAN. Additionally, we will plot 
the phase diagrams for each method, as well as the discrepancy between the 
prediction and the VISUAL NASTRAN simulation results. This will better validate 
the precision of the Matlab code. 
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Figure 50. 
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Figure 53. VN 4-D - Phase diagram - a=TT/18 - 6o=Tr/10 



Figure 54. MATLAB - Phase diagram - a=TT/1 8 - 6o=Tr/1 0 
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Figure 62. 
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Figure 63. VISUAL NASTRAN - a=TT/6 - 0 o=tt/1 0 
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Figure 64. MATLAB simulation - a=TT/6 - 0o=Tr/1O 
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Figure 66. 
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Figure 67. MATLAB simulation - a=n/6 - 6 o=tt/10 
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Figure 70. Phase diagram 
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Evaluating the above set of figures, we conclude the following: 

• The variable frequency of oscillation, as well as when 
the string switch takes place, is relatively well predicted. The error between the 
predicted and actual performance causes spikes in the angular velocity 
discrepancy plots. 

• In terms of the prediction of the angular displacement 
and velocity of the oscillation, the discrepancies are caused, on the one hand, by 
the failure to precisely predict both the location in time of the switch points and 
the energy loss at every one of them. On the other hand, it was not possible to 
obtain the exact equations used in VN 4-D to calculate these parameters. 
Therefore, a comparison with the VN 4-D software results does not provide a 
precise appreciation of the actual accuracy of the results. 

• For smaller values of characteristic angle a, which 
equals less energy loss at the switch point, the prediction discrepancy is 
minimized. 

b. String Tension Comparison 

For the string tension prediction, equation (2.46) was used. The 
Matlab code had already computed the angular velocity and displacement inputs. 
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Figure 71. Tension comparison a=TT/1 8 - 6o=Tr/10 



Figure 72. Tension comparison a=n/12 - ^o='^/^0 
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The above plots lead to the following conclusions: 

• It can be seen that, for larger values of a, which correspond 
to extensive energy loss at the switch point, VN 4-D reveals time periods when 
the tension on both strings is zeroed. This corresponds to an instantaneous 
jump. Because the mechanism of energy restitution of the ropes used in VN 4-D 
is not known and that there is the possibility that the previously assumed zero 
parallel velocity component is not completely eliminated, the resulting jump after 
impact could be expected. 

• The tension predicted by VN 4-D shows spikes of irregular 
magnitude at the switch points. The existence of these features can be related to 
the previously presented explanation. Although they do not show a pattern, the 
fact that a previously slack string becomes taught within an infinitesimally small 
time duration can provoke this peak tension magnitude. This corresponds to the 
impulse, which zeroes the parallel velocity component. 
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• The Matlab code tension prediction coincides with the VN 4- 
D only at the instances when the mass velocity is zero (indicated with a purple 
circle in Figure 72). Therefore, the VN 4-D tension prediction is not in accordance 
with equation (2.46). To understand this unexpected result, a model of two 
standard pendula -- one with a string and the other with a massless rod -- was 
used to calculate the tension during free oscillation with an arbitrary initial 
displacement. The results revealed that the standard rod pendulum’s tension 
agrees with equation (2.46), while the standard string pendulum’s tension does 
not. Moreover, an effort was made to understand the pattern that provided the 
string tension results. After some trial and error efforts, the equation that provided 
an almost identical tension prediction was 


_ -mOH , ^ 

T = —-—+m^cos6' 


(2.57) 


Efforts were made to understand the above equation, but, unfortunately, 
the MSC Software no longer supported the specific software. Consequently, 
further tension comparisons with VN 4-D will not be performed and, also, the 
Matlab results will be presented without validation from the VN 4-D results. 


3. The Real 2-Pendulum with a>45 

In the case that a=45°, when applying equation (2.55), we realize that for 
r = 0 all the kinetic energy of the pendulum would be absorbed by the string 
during impact. Thus, the oscillation would be terminated. The above can be 
visualized from the following graphs, where VISUAL NASTRAN and MATLAB 
simulations are used to provide the response of the 2-string pendulum with 
0=45°. 
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Figure 75. MATLAB simulation - a=TT/4 - 6o=Tr/10 


68 






































9' (rad/sec) 


0.3 


0.2 

0.1 

® 0 
to 

- 0.1 

- 0.2 

-0.3 

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 



t (sec) 

Figure 76. Prediction discrepancy 


2 

1.5 

1 

0.5 

0 

-0.5 

-1 

-1.5 

-2 

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 

t (sec) 

Figure 77. VISUAL NASTRAN - a=TT/4 - 0 o=tt/1 0 
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Figure 78. MATLAB simulation - a=TT/4 - 0o=Tr/1O 


1.5 


0.5 


CD 

&0 


-0.5 

-1 

-1.5 

-2 


0 0.05 0.1 0.15 


0.2 0.25 0.3 

t (sec) 


0.35 0.4 0.45 


Figure 79. 


Prediction discrepancy 
70 















































0'(racl/sec) 0' (rad/sec) 


z 



T- I 4 =I 

1 1 1 1 

1 1 1 1 


1.5 



--1 

1 1 

1 1 

-h-1- 

1 1 

1 1 


1 

— 


--1 

1 1 

1 1 

-h-1- 

1 1 

— 

0.5 

n 

— 


-- -1 

1 1 

1 1 

-h-1- 

1 1 

1 1 

1 1 

1 1 

— 

0.5 

— 


_^_ 

1 1 

1 1 

— ; — 

— 

-1 

— 


r 1 

1 1 

1 1 ^ 

-r -I 

1 1 

— 

1.5 

— 


r 1 

T -i 

1 1 

1 1 

— 

-2 


_ 

1 1 1 1 

1 1 1 1 

E_I_I_3 



- 1 - 1 -T- 1 - 1 - r ' 

-0.3 -0.2 -0.1 0 0.1 0.2 0.3 

0 (rad) 


Figure 80. VISUAL NASTRAN - a=n/4 - eo=n/10 
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Figure 82. Tensions a=TT/4 

Although the limiting case of a=45° can be described by equations (2.18) 
and (2.55) (assuming only that r = 0 ) if a>45°, the motion of the mass will not be 
governed by linear or nonlinear equations of motion, but only by its interaction 
with the strings. 



Force and velocity diagram for a>45° 
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Figure 83. 











































According to Figure 83, the linear velocity before the switch point Vb (red 

vector) is composed of two components, Vbi and Vbi (turquoise and green 
vectors, respectively). At the instance when the right string becomes tight, the 

tension of the string (blue vector) exerts an impulse affecting Vb 2 . Irrespective to 

the value of the restitution coefficient that governs this impact, Vb\ cannot 
promote the circular motion of the mass with respect to the right string due to its 
direction. Contrarily, the left string, which becomes slack after the switch point for 
a<45°, continues to be tight and exerts, along with the right string, an impulse on 

the mass opposing Vbi (orange vectors). Thus, either the mass will bounce or it 
will come to rest. This depends on the restitution coefficient. Consequently, if we 
preserve the assumption that r = 0 for all mass-string impact interactions, after 
the mass reaches the switch point, it comes to rest. 

E. THE DRIVEN PENDULUM 

Extending the analysis performed in the previous sections, we should 
include the case of the pendulum in which the supporting point, or points, 
undergoes motion due to forcing. The moving support induces complexity to the 
response of the pendulum that may not end in oscillatory behavior [6]. 

1. The Standard Driven Pendulum 

Consider a standard pendulum which, at an arbitrary moment of the 

TT 71 

oscillation, has an angular displacement 6, where -—<9< — . We define a 

world system of axis where x-axis points to the right and y-axis points upwards. 
An arbitrary motion drives the support point of the pendulum, whose acceleration 
at that moment has both a horizontal and a vertical component. 
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Figure 84. The standard driven pendulum 


According to the previous analysis, the equation of motion for the undriven 
standard pendulum is either equation (2.6) or (2.2). The effect of the horizontal 
component of the support acceleration is an inertial force on the mass 

Fx = -mX . Similarly, for the vertical component, the inertial force is Fy = -mY. 
We notice that Fy has the same direction and orientation as gravity. These two 

forces contribute to restoring force with two additional components depending on 
angle 6. Therefore, the new restoring force is 

Fr = (mg+ mY)sm0-mX cos 0 (2.58) 

Therefore, the equation of motion becomes 

(9 + ^^^^^ sin6>-^cos^ = Q (2.59) 


74 






Finally, from equation (2.59) we can draw the conclusion that, due to the 
horizontal acceleration, the pendulum will initiate oscillations even without the 
existence of initial displacement or initial velocity. Thus, for 6=6o=0, equation, 

X 

(2.59) becomes 9 = — . Finally, the above equation is valid for the range of 6 

already mentioned with one exception; when a jump occurs. The jump behavior 
of the driven pendulum will be analyzed in a subsequent section. 

2. The Driven 2-Pendulum 

Based on the previous analysis and equation (14), we can derive the 

TT 71 

equation of motion for the driven 2-pendulum as follows for <9-^a< — \ 


9+ sin(6> + a)-^cos(6>-}-a) = Q 


9 + (sin 9 cos a + sign9 cos 9 sin ~ ^ (sign9 cos 9 cos a-sin 9 sin a) = 0 


9 + ( ^^^^^ cos a + ^sin a) sin 9 + ( ^^^^^ sin a -^cos a)sign9 cos 6^ = 0 (2.60) 
Again, the validity of equation (2.60) is questioned when a jump occurs. 
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Figure 85. The driven 2-pendulum 

a. String Tension investigation of the Driven 2-Penduium 

Assuming that the support of the 2-pendulum is accelerated in the 
direction shown in Figure (85), the mass will leave the equilibrium position only if 
72=0 N. For this to occur, the vertical component of 7f should be able to 

withstand the virtual weight of the mass, which is Wv = m(Y -h g). We can easily 

prove that the magnitude of the string tension about which the mass will oscillate 
is 

T = m{0H + X?,m{9+a) + Y co^{9 + a) +gco^{0 + a)) (2.61) 

TT 71 

which is valid for all -—+a<9<—-a. Nevertheless, close attention should 

2 2 

be given in terms of the sign of angles B and a If inconsistencies are to be 

avoided, angle a should always have the same sign as angle 6. A quick look at 
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Figure 86 shows us that a positive X would increase the tension on the string 
while a negative X would decrease it. Nevertheless, if the mass were oscillating 

about the other string, a negative X would increase the string tension, and so 
on. The above facts are predicted by equation (2.61). 



Figure 86. 2-pendulum’s force and acceleration diagram during oscillation 

Nevertheless, equation (2.61) has one limitation: the tension can 
only take positive values. When the string tension is zeroed, a jump occurs. That 
behavior will be analyzed in a subsequent section. 

In addition, it would be useful to know the magnitude of the tension 
exerted on the strings while the mass remains in its equilibrium position. While 
the mass is in that position, the sum of the vertical components of Tf and 72 
should be equal to the virtual weight of the mass Wv, while the difference of their 

horizontal components should be equal to mX 

77 





co^a(Tl + T2) = m(Y + g) and ^ma(Tl-T2) = mX 

solving the above system of two equations 

7^1 _ ^ ^ ^ ^ ^ 

2 ^ cos a sina"^ 

j2_m^ (Y + g) X ^ 

2 ^ cos a sina"^ 


(2.62) 


b. Stability of the 2-Penduium 

While equation (2.60) describes the 2-pendulum’s oscillation in the 
absence of non-zero initial conditions, the oscillation will not be initiated under 
any value of horizontal acceleration. 


0 = - 


Manipulating equation (2.60) with zero initial conditions, results in 
^^^^^sincr-h^coscr. The vertical components of the restoring force 


tend to drive the mass to move downwards. This motion is constrained in this 
situation since the mass is supported by both the taught string (an exception 

exists for the case that Y is negative and ^-hF<0). Therefore, an angular 


acceleration will be induced only if the horizontal component manages to 
overcome the condition and 


is fulfilled. 


X 

tana 


Y>g 


(2.63) 


Another approach, which would lead us to the same conclusion, is 
analyzing the string tension. Assuming that the support of the 2-pendulum is 
accelerated in the direction shown in Figure 86, the mass will leave the 
equilibrium position only if 72=0 N. For this to occur, the vertical component of 7f 
should be able to withstand the virtual weight of the mass, which is 
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Wv = m{Y + g). In the case of the 2-pendulum, which is about to leave its 

equilibrium position 6 = 0 and from (2.61), for the mass to be lifted up we must 
have 

T cos a = m(X sin a -f 7 cos a + g cos a) cos a > m{Y -t- 
X sin a cos a> g sin^ a + Ysin^a 


tana 

Finally, from equation (2.63), we can draw the conclusion that 
choosing larger values for a would make the 2-pendulum more stable. 
Nevertheless, this advantage has a considerable drawback: the increase in string 
tension associated with larger characteristic angles. In Figure 87, one can 
observe both the limiting value of horizontal acceleration, which allows the 2- 
pendulum equilibrium, and the maximum tension associated with this horizontal 
acceleration for a range of characteristic angles (1°<a<80°). 
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Figure 87. Horizontal acceleration inducing instability versus associated maximum 

tension 
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From the above figure, it can be seen that, although the necessary 

X would break the equilibrium, which increases steadily for the first half of the 
range of a, the related tension increases in a slower pattern. This difference in 
increasing pattern is caused by the analytical equations that describe them. 
These equations are 


^i.m = (^ + 4m)tana 

T 


_ +^max) I 


( 2 . 64 ) 


cos a sin a 

One method to evaluate the optimum choice of a would be to define 
a ratio of acceleration over the maximum attained tension. According to the 
previous equations, this ratio would be R = sina. 


Nevertheless, this method is not effective. The values of the 

required X to disturb equilibrium are very large for the second half of the a 
range, and even for some values of the first half. A more realistic approach is to 
choose a value of a (for example 20°) as a reference value, Qref, and, then, 
compare the maximum tension, Tmaxref, with the one exerted on the strings of the 

2-pendulums of larger a, but still associated with the same X^^j-. 
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Figure 88. a=20° reference angle 


The above plots reveal that, in a choice of a greater, but less than 
52°, a would be a better choice. 

Another approach is to find the angle a which will result in the 
minimum exerted tension for a combination of vertical and horizontal 
accelerations. In the following figures, we can visualize the variation of optimum 

a when X and Y vary from zero to g along with associated tension. The results 
are validated with the constraint that the 2-pendulum should be stable for the 

combination of a, X , and Y. 





I 

I 

I 

I 

' ^ 

I 

I 

-*-^- 

-h- 

I 

I 

I 

I 

I 

I 


____J_. 

I 

I 

I 

I 

_^__ 


81 




























Tension (N) 0 a (degrees) 


Characteristic angle giving the minimum tension 



89. Optimum characteristic angle investigation for combinations of 

accelerations 


Minimum tension associated with optimum characteristic angie 



Figure 90. Tension related to optimum a 
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A first step to take advantage of the above plots and analysis would 

be to choose the maximum expected X and Y for our design and, then, iterate 
for the optimum a in terms of minimum tension. The second step would be to 

verify that the design would be stable for a new range of X and Y composed of 
negative values because Y tends to destabilize the 2-pendulum. In the following 
example, the range of -0.3^ <7 <0.7^ is chosen to challenge the stability of 
the selection of a=38° designed to provide the minimum tension for the 
X = g,Y = g combination of maximum accelerations. 

As depicted in the following figure, there exists a range of X and 
Y, where the driving design constraint is not the desired minimum tension, but, 
rather, the stability for the pendulum. Therefore, for X = Y = g , \Ne should select 
another characteristic angle according to the stability criterion. 


Design Driving Constraint 



Figure 91. 


Design driving factor constraint versus objective function minimization 
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Characteristic angie for combined constraints 



Figure 92. Characteristic angle selection for combined constraint and objective 

function minimization 

Therefore, from Figure 92, we can decide which should be angle a. 
This would be based on the worst expected negative vertical acceleration. For 

example, a combination oi X = g,Y = -0.3g would require a=55°. 

In addition to the above manipulations, another design method 
would be the use of basic optimization techniques. In terms of optimization, the 
problem is a minimization problem, i.e., minimizing the maximum exerted tension 

(^max = ) with one objective variable (characteristic angle o) 

cos a sm a 

and subject to boundary conditions {a>0,a> arctan(-^—) ). A Matlab code 

g + Y . 

o min 

was developed to yield the optimum characteristic angle for given maximum and 
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minimum accelerations. In the following figure, we can visualize the results for an 
arbitrary input combination. 


inaxXddot=0.4g - inaxYddot=0.3g - ininYddot=-0.3g 



Figure 93. Optimum selection for arbitrary design inputs 

c. Jumps of the Driven 2-Pendulum 

The tension magnitude provided by equation (2.61) denotes that, 
under certain conditions, string tension may become zero. Keeping the 

TT 71 

assumption that -—-\-a<9<—-a , the only terms that can take tension off the 

string are X^mO and Ycos9. Therefore, the restriction that ||^+a||>y for a 

jump to occur, which is valid for the unforced 2-pendulum, is unnecessary. To 
predict the motion of the mass after the jump, we should follow a certain 
procedure. 
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Jump 

Jump 


end 

initiation 




Figure 94. Jump of the 2-pendulum 


Firstly, just after the jump, we should track the position with respect 
to a global coordinate system, xo and yo, and to the linear velocity magnitude, u, 
of the mass. With the above initial conditions, the mass will follow a projectile 
path. Thus, the equations of motion are 

X = Xq+UoICOS^ 

. ^ 1 2 (2.65) 

y = yQ + UQtsm^--gt 

Secondly, after the jump, we should track the positions, velocities, 
and accelerations of the two support points. Depending on the type of motion 
applied to the support, we should especially track the position of the support 
points during the jump for each time step. 

Finally, we should track the distance between the mass and each 
one of the support points. When this distance becomes equal to the string length. 
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it is dynamically possible for the oscillation to proceed again. The methodology 
used to predict the oscillation characteristics after the end of the jump is similar to 
the one used while analyzing the pendulum’s behavior at the switch points. 
Firstly, the jump takes place while both the strings lengths are less than /. 
Consequently, the first step is to identify which string becomes taught first. This 
determines if angular displacement has occurred after the jump. Secondly, the 
magnitude of post-jump angular displacement has to be calculated. Knowing the 
mass positions just before the end of the jump, x and y, as well as each of the 


string’s support point positions, =Z-/sina and = Z-h/ sin a, the 
total angle defined by the taught string and the vertical, cp is 


(f) = tan( 
(f) = tan( 


^ right ^ 


v 


right 


y 


^), 6>0 


z-Z 


left 


^14.-y 


),d<0 


( 2 . 66 ) 


Therefore, 


0 = (p-a,6>0 
6 = -^ + a,d <0 


(2.67) 


Moreover, at the moment when one of the strings becomes taught 
one faces again a mass to string impact. Knowing the linear velocities of the 
mass just before the impact. 


=MqCOS^ 

Uy =Uq sin^-gt 


( 2 . 68 ) 


and the horizontal - vertical velocities of the support, X and Y, the mass will 
undergo impact with the taught string. This is governed by the same restitution 
coefficient assumption that was explained in Chapter II, Section D, number 2. In 
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this case, the impact velocity components are the relative velocities between the 
mass and the support, =u^-X,Uyi,^p =Uy-Y. Consequently, the following 
two equations are used to obtain the post-jump angular velocity 


^ = \ cos (j) + sin ^), 6' > 0 

<^ = j i-Uxi^p cos ^ - Uy^,^p sin ^), < 0 


(2.69) 


Finally, knowing the post-jump angular displacement and velocity 
initial allows the continuation of the oscillation subject to the new initial 
conditions. 


3. Evaluation of Results 

In the following two sets of figures (Figures 95-99 and Figures 100-104), 
two examples of the response prediction for forced oscillation with jumps are 
presented. A Matlab code Program_3.m was developed to predict the position of 
the mass while the supports are subjected to horizontal, for the first set, and to 
horizontal and vertical motion, for the second set, under conditions that will force 
a jump to occur. The motion that was applied to the support midpoint can be 
visualized in Figures 95 and 100. It can be seen that in order to produce a jump, 
i.e., zero tension on both strings, we applied a sudden stop to a motion that was 
previously under the effect of variable acceleration. This violent stop was enough 
to produce a distinguishable jump both in VN 4-D and in Matlab. The blue solid 
lines represent the Matlab predictions for the mass positions and angular 
displacement, while the red dashed lines represent the VN 4-D predictions. 
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Figure 95. Support midpoint trajectory 
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Figure 96. Mass trajectory 
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Figure 97. Angular displacement 


1 

0.5 
0 

-0.5 

-1 

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

t (sec) 

1 

0.5 
0 

-0.5 

-1 

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

t (sec) 

Figure 98. Mass absolute position 
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Figure 99. 


Tension variation 
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Figure 100. Support midpoint trajectory 
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Figure 101. Mass trajectory 



Figure 102. Angular displacement 
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Figure 103. Mass absolute position 



Figure 104. Tension variation 
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On the one hand, the positions of the support midpoint (Figures 95 
and 100) and the total mass positions (Figure 96 and 101) are calculated to 
coordinate the axis system with its origin positioned on the support midpoint Xg- 
Yg. On the other hand, the absolute mass position is calculated wrt a coordinate 
axis system with its origin at the equilibrium position of the mass following the 
support motion (Figures 98 and 103). Moreover, the angular displacement is 
presented where, during the jumps, the value is set to zero (Figures 97 and 102). 
Finally, as explained in Chapter II, Section D, number 2.d, the strings’ tensions 
are plotted without including the VN 4-D results. 

The discrepancy that can be observed between the Matlab and VN 
4-D results does not seem to be related to the complexity of the induced support 
motion. After a more extensive study of the VN 4-D results (tabulated data 
exported from the software interface), it can be seen that, in some of the jumps, a 
nonlinear change in the mass velocity exists at the exact moment of the jump 
initiation. Introducing this exact nonlinearity in the Matlab code yields the same 
results with VN 4-D. 

4. Convergence Evaluation 

An essential and final step that must be performed before ending the 
analysis of the driven 2-pendulum is to examine the convergence of the results 
for different choices of time steps. From the initial analysis of the free 2- 
pendulum to the simulation of the driven 2-pendulum with jumps, the integration 
time step selected for the numerical solution of the nonlinear equations of motion 
was df=0.001 sec. Therefore, a basic evaluation procedure of the examples will 
be presented to verify that the selected time step offers convergence of results 
over other choices of larger time steps. 
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Figure 105. Convergence evaluation for forced oscillation - no jumps 



Figure 106. Convergence evaluation for forced oscillation - jumps 
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In Figure 105, the convergence is tested with both vertical and horizontal 
support motion. A case with jumps is presented in Figure 106. From both figures, 
one can easily conclude that gradually increasing the time step dt will provoke a 
solution divergence. While the discrepancy for df=0.005 sec or df=0.01 sec (up to 
ten times larger than the initial choice) will cause the solution to change slightly, a 
further increase in the time step size may yield solutions with larger 
discrepancies than the original. Especially in the jump cases, it can be seen that 
some jumps are incorrectly predicted (in terms of time of jump initiation and 
termination) or even skipped if the time step is large enough. 
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III. THE FOUR-STRING PENDULUM 


A. THE 4-PENDULUM DESIGN ANALYSIS 


To extend the previous analysis of the 2-pendulum to more complex 
designs, the performance of a four-string pendulum will be analyzed. The four 
strings have the same length and their supports are symmetrically placed wrt the 
mass. All the strings are attached to the center of mass. The following figures 
represent the model that was used in VN 4-D and offers four views of the design: 



Figure 108. 4-pendulum bottom and side view 
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Such a design is one of the more representative simplifications of the real 
world suspension designs. Multi-rope suspension designs differ from each other 
on the number of ropes that are used, the positioning of their support points, and 
the selection of their attachment positions on the suspended load. 

The purpose of this section’s analysis will not be to investigate the 
swinging attitude of such a pendulum, but, rather, to study its ability not to swing 
when accelerating motion is applied to its support. Moreover, the string tensions 
will be analyzed to try to draw the most accurate picture of the design 
characteristics that should be considered to achieve a desirable 4-pendulum 
design. 



Figure 109. 4-pendulum design characteristics 


In the above figure, the black sphere represents a Ikgr mass and the 
black lines are the strings of the pendulum. The global coordinate system of axis 

is presented. It dictates that the positive X-axis represents the forward motion; 
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the positive Z-axis represents the side motion towards the right; and the positive 
Y-axis represents the upward motion of the support. Since we are now dealing 
with a 3-D analysis, instead of 2-D for the 2-pendulum, there are two distinct 
characteristic angles that have to be defined. In Figure 109, the angle defined by 
the two green lines will be called the forward principal characteristic angle, cpi, 
while the angle defined by the two turquoise lines is the side principal 
characteristic angle, (p 2 . 

1. String Tension Investigation of the Driven 4*Penduium 

The 4-pendulum is assumed to be subjected to a forward, lateral, and 

vertical support motion characterized by X, Z, and Y accelerations, 
respectively. To formulate the tension equations for each individual string, the 
superposition rule will be used. Let us name the strings: string 1, the one closest 
to the global system of axis; and strings 2-4, counting clockwise looking the 
pendulum from the top (Figure 109). 
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Figure 110. 4-pendulum front view tension breakdown 


Initially, it will be assumed that X = 0. Thus, the analysis will be similar to 
the one completed for the 2-pendulum. In the figure above, the side view of the 
4-pendulum shows that the strings, represented as dashed turquoise lines, lay on 
the same plane with the turquoise non vertical lines of Figure 109. If one of these 
dashed turquoise lines substituted each side pair of strings, the principal tension 
Tp in each of those will be in accordance with equations(2.62). 
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(3.1) 


Tp2\ 

Tpii 


m^ {Y + g) ^ Z ^ 

2 ^ cos ^2 sin ^2 


m^ {Y + g) 

2 cos ^2 


Z 

sin ^2 


) 



Looking the same model from the left-side view, Tpii can be analyzed 
into two components, Tp 2 \\ and Tp2i2, since they correspond to string 1 and 2, 

respectively. Let us name the angle that each of those tension components forms 
with the resultant tension, (Ppartiai 2 - These two components are equal and their 
magnitude is 


Tp2\l = Tp2\2 


Tp2\ 

'2 cos ^ partial! 


(3.2) 
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Then, combining equations (3.1) and (3.2), 

m 


Tp2\l = Tp2\2 = 




-) 


^08^2 Sill ^2 

Performing the same analysis for strings 3 and 4 we yield to 


(3.3) 


7/9223 = Tp22A = 


m 


-( 


(5^ + ^) 


-) 


(3.4) 


cos^ 2 sin^2 

If we rename the tensions using simpler symbology to T1, T2, T3, and T4, 

respectively, and perform the same methodology for the case that Z =0, the 
yielding results are 

m 


T\ = TA = 


4cos^ 


T2 = T3 = 


partiall 

m 


Aco^(/) 


partial 1 


cos^j sin^j 

^ (y+8) X ^ 

^ cos^4 sin^j 


(3.5) 


Introducing three new and more generalized variables H, as the vertical 
position of the mass measured along the negative Y-axis, Lfront, as the frontal 
distance between the support points, and Lside, as the side view distance 
between the support points, the principal characteristic angles can be expressed 
as follows: 


/ , Ls idc X 

^i=arctan(— 

, (3.e 

^2 = arctan(^^^^^) 

Additionally, let us define Lpiand Lp 2 , the lengths of the pseudo-strings 


of the principal tensions (green and turquoise non-vertical lines in Figure 101), 
and /, the length of each string. So, 
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(3.7) 



=ar':cos(^) 


It can easily be proved, using the two previous equations, that 

^partial! A ^ partial! COS ^2 



Therefore, combining (3.3), (3.4), (3.5) and (3.8), the tensions exerted on 
the four strings for every combination of forward, lateral and vertical 
accelerations are 

j^_m^ (Y + g) ^ 2X ^ 2Z ^ \^2 I ^ Lfront ^ ^ ^ Lside ^ 

4 H Lside Lfront \ 2 2 

^2-^/ (^ + g) 2X ^ 2Z 1^2 I r Lfront .2 , r Lside ^ 

4 H Lside Lfront v 2 2 

(3.9) 


^^_m^ {Y + g) ^ 2X 2Z ^ 
4 H Lside Lfront 



^^_m^ {Y + g) 2X 2Z 
4 H Lside Lfront 




2 , .Lfront.2 , rLside.2 


-r+(- 
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This set of equations can be used, as already done for the corresponding 
one of the 2-pendulum, to predict the stability of the 4-pendulum, i.e., at which 
point that at least one of the string tensions will be equal to zero. 

2. Tension Analysis Validation 

Aiming to validate the above set of equations, the VN 4-D model was 
used. The setup of the model was complicated due to the required level of 
dimensional precision. Finally, the method that provided the most reasonable 
results (absence of nonlinearities in tension) was to combine two separate runs 
of the software - one with no forward acceleration and one with no lateral 
acceleration - and then superimpose the results. In this specific model m=1kgr, 
the frontal distance between the support points is /./ronf=2m; the side distance is 
Lside=^^f]■, and the vertical position of the mass is H= 2 m. Finally, the applied 

accelerations are X =2f(m/sec^), Z = 4f(m/sec^), and Z = 0 (t in seconds). 



L_^^^_I_ I _^_ 

0 0.1 0.2 0.3 0.4 0.5 0.6 

t (sec) 


Figure 112. Strings 1 and 2 tension validation 
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Figure 113. Strings 3 and 4 tension validation 


B. 4-PENDULUM DESIGN OPTIMIZATION 

The procedure of design optimization for the 4-pendulum requires that one 
defines the objective function, the design variables, and the constraints - just like 
it was done for the 2-pendulum. Once again, the objective function is to minimize 
the tension that will be exerted on a string for a given set of accelerations. The 
design variables are the distances H, Lfront and Lside. 

According to (3.9) set of equations, the Lfront, Lside, and H substitute 
and ^2 the design variables. Nevertheless, we can constrain the value of H 

with an equality constraint, instead of an inequality one, to be able to obtain a 2- 
D visualization of the objective function. 

In the following example, the known values are /-/=1m, 
lfront ■ = Lside ■ = 0.2m, and 
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Figure 114. Maximum exerted tension for arbitrary accelerations 


106 



















Minimum Tension (N) 


O) 

c 

a; 


d) 

-q 

'(/) 


— 

11 

1 





I 

I 

I 

I 

I 

I 

I 



7 .? 


I 

00 

T- 

- 1 

)- 



I 

I 

I 

I 

I 

I 

% ■ / 

S' / 1 

o ' 



N) 

I 

I 

I 

I 

I 

I 

1 


o 

I 


I 

I 

I 

I 

I 

I 

I 

1 I 

V 1 _ AO 

7 ^ 1 A2 

I 

14 

16 A8 -2 

aa 

0 - -22^=^ 

A6 

>4= 

I 

0 1 2 3 4 5 6 


Frontal length (m) 


Figure 115. Maximum tensions exerted for arbitrary accelerations 



Figure 116. Visualization of designs that fulfill the stability constraint 
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Figure 117. Final visualization of maximum tensions for arbitrary accelerations 

In Figure 117, we can visualize the acceptable Lfront and Lside 
combinations and at the same time find the optimum solution to our problem. 
Clearly the optimum solution is Lfront= Ls/de=2m. This solution can also be 
verified using the fmincon optimization Matlab function, which also yields the 
same result. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


A. SUMMARY AND CONCLUSIONS 

The following table summarizes the basic equations developed throughout 


this thesis. 


A/A 

Description 

Equation 

Equation 

No 

1 

Free 

oscillation 

of 

2- 

pendulum 

^ + y sin 6cos a + y cos 0 sin a sign 9 = Q 

subject to 

0^^ = 0_,^ism{9O-2a)f + (r cos(90 - 2a)f 

(2.18) 

and 

(2.55) 

2 

SCPA 

method 

prediction 

f 2P ^ 

6' = Acos CD^ + —-— t 

? s 

cOq =yCosa 

P = ysina 

(2.32) 

and 

(2.24) 

3 

Energy 

method 

prediction 

^ + Fysin^ = 0 

F =\^^mgsm[a + 0)d0 l\Q^mgsm0d0 

(2.39) 

and 

(2.38) 

4 

Driven 

2- 

pendulum 

X Ag+y) X . ... 

+ —^cosa +—sma)sin6' + 

^ sin a - ^cos a)sign0cos 0 = 0 

(2.60) 
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stability 
criterion for 

driven 

2 - 

pendulum 




tana 




(2.63) 


String 
tensions 
for driven 

2 - 

pendulum 
while in 
equilibrium 
position 


j^_m^ {Y + g) ^ X ^ 
2 ^ cos a sina"^ 
^^_m( (7 + g) X ^ 
2 cos a sin a 


(2.62) 


String 
tensions 
for driven 

4- 

pendulum 
while in 
equilibrium 
position 


_mAY + g) 2X 2Z \ 2 ^(Lfront. 2 ,Lside. 2 


zSs" 


^^-^( (F + g) 2X ^ 2Z ^ ^2 I ^ Lfront ^2 , ^ Lside ^2 

4 H Lside Lfront \ 2 2 


_mAY + g) 2X 2Z \ 2 ,(Lfront. 2 , (Lside. 2 




4 H Lside Lfront 




i \ '} / i-diy n^c. N •) 

-) +( - ) 

2 ^ ^ 2 ^ 


_mAY + g) 2X 2Z 2 ^(Lfront. 2 ^ (Lside .^2 




4 H Lside Lfront 




2 ^ 2 ^ 


(3.9) 


Table 2. Equations summary 


Based on the analysis performed throughout this thesis, including both the 
analytical formulation of the mathematical equations and their validation 
performed with the aid of VN 4-D, we can draw the following conclusions: 
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• The 2-pendulum’s configuration of two strings per mass 
introduces two distinct nonlinear oscillation characteristics. The first one results 
from the discontinuity in its stiffness around the equilibrium position, i.e. an 
instantaneous jump in the value of the restoring moment around zero. Four 
different methods, were used to incorporate this discontinuity in analytical 
equation formulation: Numerical solution with ode4 5, exact solution for a time 
increment, SCPA method and Energy method. 

• The second nonlinear oscillation characteristic is related to 
the energy dissipation that takes place at the switch point. The dynamics that 
govern this behavior are correlated to the dynamics of a standard collision of a 
mass on a rigid wall, with specific velocity direction constraints before and after 
the collision. The idealized model is based on an instantaneous impact and total 
energy loss related to the mass’ velocity component parallel to the taught string’s 
direction. This assumption is validated by the VN 4-D simulation. 

• Another characteristic that was analyzed was the stability of 
the 2-string pendulum in forced oscillation. This string configuration results in the 
absence of oscillatory behavior if the mass is initially at the equilibrium position 
and the driving horizontal and vertical accelerations do not exceed values 
determined exclusively by the characteristic angle a. This observation allows the 
rational design of support systems under given support motion characteristics. 

• The general behavior of the driven 2-pendulum was 
investigated including the case of induced jumps. Specific examples 
demonstrated the highly non-nonlinear behavior of the 2-pendulum when 
extreme acceleration values of the support cause the strings to go slack. 

• Finally, the analysis of the 2-pendulum was extended on the 
4-pendulum configuration, which is a more common support design in practice. A 
typical optimization study was presented in order to visualize the practical design 
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advantages that such an analysis offers. It was shown that it is possible to 
reduce both motions and cable tensions by appropriate selections of string 
angles. 

B. RECOMMENDATIONS 

Further development of the analysis presented in this work could be in the 
application of angular support motion, both for the case of forced oscillations of 
the 2-pendulum and the trade off analysis. This work would provide a more 
realistic representation of the support motion analogous to the motion of the ship. 
Further extensions to this work by considering the effects of string elasticity are 
also recommended. 
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APPENDIX: 


COMPUTER PROGRAMS 


%Program_l.m 


"6 

"6 

%Computing the response of 
methods: Numerical solution 
"6 
"6 


free 2-pendulum using three distinct 
, SCPA method and Energy method 


clear all 
clc 

close all 

global a omega2 A omega P maxyl Thetao t increased_restoring 
THETAo=pi/10 %initial angular displacement 

THETAo_dot=0 %initial angular velocity 


a=input( 'Give the characteristic angle a ' ); %2-pendulum's 
characteristic 


g=9.81; 
11 = 1 ; 


l=ll/cos(a) 
omega2=g/l; 

omega2o=omega2*cos(a) ; 
tf=3; 


%angle 

%acceleration of gravity 

%equivalent string length which equals to the 
%simple pendulum string length or the 
%perpendicular distance of the mass from 
%the plane of suspension 
%string length of complex pendulum 
%square of the natural frequency of the 
% standard pendulum 

%square of the pseudo natural frequency of the 

%2-pendulum 

%final time 


"6 

"6 

%Numerical solution of the differential equation of motion 

"6 

"6 

options = odeset( 'RelTol' ,le-6, 'AbsTol' , [le-6 le-6]); 

[t,x]=ode45 ('Ffull', [0:0.01:tf],[THETAo,THETAo_dot],options); 
[t,x(:,l) ] ; 

Xl=x; 


figure (1) 


hold on 

title( '\alpha=\pi/6 ' ) 

plot(t,X(:,1), 'g' ) 
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maxyl=max(x(: , 1)) ; 
maxy2=max(x(:,2)) ; 


"6 

"6 

%SCPA method 

"6 

"6 


A=THETAo; 

B=sqrt(2*g*(11-cos(THETAo+a))/I); 

t=0:0.01:tf; 

[xap,yap,Tapprox,Uo,Texact,t1,Ain]=SCPA(t,a,omega2,A,omega,P,maxyl); 


plot(t,xap, '—r' ) 


Xlap=[]; 

Ylap=[]; 

for 1=1:Tapprox*100; 
t=(1-1)*0.01; 

[xap, yap, Tapprox, Uo, Texact, 11,Ain]=SCPA(t,a,omega2,A,omega,P,maxyl); 

Xlap(1)=xap; 

Ylap(1)=yap; 
end 

"6 

"6 

%Energy method 

"6 

"6 

syms X Theta phi increasing_factor square_nat_freq phi t Initial_theta 

Initial_theta_dot 

y=int('sin(x) ',x,0,Theta); 

z=int( 'sin(phi+x)' ,x,0,Theta); 

z=subs(z,phi,a); 

z=subs(z,Theta, maxyl) ; 

y=subs(y,Theta,maxyl) ; 

increased_restoring=z/y; 


[t,y]=ode45 ('E', [0:0.01:tf],[THETAo,THETAo_dot],options); 
[t,y(:,l) ] ; 
plot(t,y(:,1), '-.b ' ) 

Yl=y; 
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t=0:0.01:tf; 


method' ) 


hold off 


xlabel ( 't (sec) ' ) 
ylabel( '\theta (rad)' ) 
grid on 
box on 

legend (' ODE45 'SCPA method', 'Energy 
axis([0 tf -1.25*maxyl 0.6]) 


Tapproxl=Tapprox; % Approximate Period of the Slowly Changing method 
Texactl=Texact; %Period from the analytical solution calculating 

% for the time period when the sign of theta is 
%constant 


figure (2) 

title( '\alpha=\pi/6' ) 
hold on 

plot(XI(1:Texact1*100,1)./THETAo,XI(1:Texact1*100,2)./sqrt(omega2o), 'g' 
,Xlap./THETAo,Ylap./sqrt(omega2o) , '—r' , 

Y1(1:Tapproxl*100, 1) ./THETAo,Y1(1:Tapproxl*100,2) ./sqrt(omega2o), '-.b' ) 
grid on 
box on 

xlabel( 'Xtheta/\thetao' );ylabel( '\theta''/\omega\o' ) 
legend (' ODE45 ' , 'SCPA method ',' Energy method') 
axis([-1.51.5-1 1]) 

hold off 
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%Subroutine Ffull.m of Program_l.m and Program_2.m 

"6 

"6 

%Numerical solution of the differential equation of motion for 
%the 2-pendulum 


function xp=Ffull(t, x) 
global a omega2 

xp=zeros(2,1); 
xp (1) =x (2) ; 

xp(2)=-sin(x(1))*omega2*cos(a)-omega2*cos(x(l))*(sin(a))*sign(x(l)); 


%Subroutine E.m of Program_la.m and Program_lb.m 

"6 

"6 

%Numerical solution of the differential equation of motion for 
%the 2-pendulum 
"6 
"6 

function yp=E(t,y) 

global a omega2 increased_restoring 

yp=zeros(2,1); 
yp (1) =y (2) ; 

yp(2)=-sin(y(1))*omega2*increased_restoring; 


%Subroutine SCPA.m of Program_l.m 

"6 

"6 

%SCPA method 

"6 

"6 

function 

[xap,yap,Tapprox,Uo,Texact,t1,Ain]=SCPA(t,a,omega2,A,omega,P,maxyl); 

global a omega2 A omega P maxyl Ain 

P=omega2*sin(a) ; 

omega =sqrt(omega2*cos(a)) ; 

tl=(1/sqrt(omega2))*acos(1/(1+(A*omega2*cos(a))/(omega2*sin(a)))) ; 
xap=A*cos((omega+2*P/(pi*omega*A))*t); 

yap=A*(omega+2*P/(pi*omega*A))*sin((omega+2*P/(pi*omega*A))*t); 
Tapprox=2*pi/(omega+2*P/(pi*omega*A)); 

Uo=-omega*(A+P/(omega.^2))*sin(omega*tl); 

Ain=-Uo/omega - (P/ (omega . ''2) ) *sin (omega*11) ; 

Texact=(4/omega)*acos(1/(l+A*omega*omega/P)) ; 
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%Program_2.m 


%Computing the response of 
methods: Numerical solution 


free 2-pendulum using three distinct 
, SCPA method and Energy method 


clear all 
clc 

close all 

global a omega2 A omega P maxyl Thetao t increased_restoring 
THETAo=0 %initial angular displacement 

THETAo_dot=l.22 %initial angular velocity 


a=input( 'Give the characteristic angle a ' ); %2-pendulum's 
characteristic angle 

g=9.81; %acceleration of gravity 

11=1; %equivalent string length which equals to the 

%simple pendulum string length or the %perpendicular distance of the 
mass from 

%the plane of suspension 

l=ll/cos(a) %string length of complex pendulum 

omega2=g/l; %square of the natural frequency of the 

% standard pendulum 

omega2o=omega2*cos(a); %square of the pseudo natural frequency of the 

%2-pendulum 

tf=3; %final time 


"6 

"6 

%Numerical solution of the differential equation of motion 

"6 

"6 


options = odeset( 'RelTol' ,le-6, 'AbsTol' ,[le-6 le-6]); 

[t, x]=ode4 5 ('Ffull', [0:0.01:tf], [THETAo,THETAo_dot],options); 
[t,x(:,l) ] ; 


%Predicting the maximum displacement based on the initial velocity 
% using energy conservation and geometry 


so=l*cos(a) ; 


117 














h=THETAo_dot*THETAo_dot*l*l/(2*g) ; 
hl=ll-h; 

THETAo2=acos(hl/1)-a; 


figure (1) 


hold on 

title( '\alpha=\pi/6' ) 

plot(t,X ( :,1), 'g' ) 


maxyl=max(x(:,1)); 


%SCPA method 


A=THETAo2; 

t=0:0.01:3; 

[xap,Tapprox,Uo,Texact,11,Ain]=SCPA_u(t,a,omega2, A, omega, P, maxyl) ; 


plot (t,xap, '—r' ) 


"6 

"6 

%Energy method 

"6 

"6 

syms X Theta phi increasing_factor square_nat_freq phi t Initial_theta 

Initial_theta_dot 

y=int('sin(x) ',x,0. Theta); 

z=int( 'sin(phi+x)' ,x,0,Theta); 

z=subs(z,phi,a); 

z=subs(z,Theta, maxyl) ; 

y=subs(y,Theta,maxyl) ; 

increased_restoring=z/y; 


[t,y]=ode45 ('E', [0:0.01:tf],[THETAo,THETAo_dot],options); 
[t,y(:,l) ] ; 
plot(t,y(:,1), '-.b ' ) 
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method' ) 


hold off 


xlabel ( 't (sec) ' ) 
ylabel( '\theta (rad)' ) 
grid on 
box on 

legend (' ODE45 'SCPA method', 'Energy 
axis([0 tf -1.25*maxyl 3*maxyl]) 


%Subroutine SCPA_U.m of Program_2.m 

"6 

"6 

%SCPA method for initial velocity 

"6 

"6 

function 

[xap,Tapprox,Uo,Texact,t1,Ain]=SCPA(t,a,omega2,A,omega,P,maxyl); 

global a omega2 A omega P maxyl Ain 

P=omega2*sin(a) ; 

omega =sqrt(omega2*cos(a)) ; 

tl=(1/sqrt(omega2))*acos(1/(1+(A*omega2*cos(a))/(omega2*sin(a)))); 
xap=A*sin((omega+2*P/(pi*omega*A))*t); 

Tapprox=2*pi/(omega+2*P/(pi*omega*A)) ; 

Uo=-omega*(A+P/(omega.^2))*sin(omega*tl); 

Ain=-Uo/omega - (P/ (omega . ■'2) ) *sin (omega*11) ; 

Texact=(4/omega)*acos(1/(l+A*omega*omega/P)) ; 
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%Program_3.m 

"6 

"6 

%Computing the response of the driven 2-pendulum for a given and 
% specific set of support motion 


"6 

"6 


"6 

9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'ZiTTT?'NTTT rM^TS- 9- 9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'2'2'2'2'2'2'2'S'2'2'2'2' 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOrl±±J_iiM±-LwiNOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 

%this m-file is valid only when the appropriate and specific data_l.m 
%and data_2.m data files exist in the same folder since those files 
% contain both the support motion info and the corresponding results 
% from Nastran for comparison and validation purposes 



clc 

clear all 
format long 

fprintf( 'make sure that you have copied the data provided by Nastran') 
datainput=input( 'if you already have press 1 ') 

if datainput==l 

fprintf (' good to go') 

else 

break 

end 


%call the data file 

dat=input (' choose the data file, 1 for latteral and vertical motion, 2 
for latteral motion') 
if dat==l 
data_l 
else 
data_2 
end 
"6 
"6 

%response before the stability is broken 

"6 

"6 


ti=0; %initial time condition 

tf=length(data(:,2))/lOOO-O.001; %final time condition 
dt=0.001; %time step 
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nt=fix((tf-ti)/dt) 

a=pi/18; 

g=9.81; 

11 = 1 ; 


l=ll/cos(a); 
omega2=g/l; 

omega2o=omega2*cos(a) ; 

r=0; 
m=l; 


%number of time steps 
%pendulation characteristic angle 
%acceleration of gravity 

%equivalent string length which equals to the 

%simple pendulum string length or the 

%perpendicular distance of the mass from the 

% plane of suspension 

%string length of complex pendulum 

%square of the natural frequency of the simple 

%pendulum 

%square of the pseudo natural frequency of the 
%2-pendulum 

%restitution coefficient 
%mass 


%initializing the acceleration, velocity and displacement vectors 


acc=zeros(1,nt+1); 
vel=zeros(l,nt + l) ; 
disp=zeros(1, nt + 1) ; 
force=zeros(1,nt + 1) ; 
Zabs=zeros(l,nt + l) ; 
Yabs=zeros(l,nt+l); 
Z=zeros(l,nt + l) ; 

Y=zeros(l,nt+l); 
Zs=zeros(l,nt + l) ; 
Ys=zeros(l,nt + l) ; 
acc_rest=zeros(1, nt + 1) ; 
%that the mass doesn't 
Zsddot=zeros(1, nt + 1) ; 
%midpoint 

Ysddot=zeros(1, nt + 1) ; 
%midpoint 

Zsdot=zeros(1, nt + 1) ; 
Ysdot=zeros(1, nt + 1) ; 
Zsdot(1,1)=0; 

Ysdot(1,1)=0; 

Zsddot(1,1)=0; 

Ysddot(1,1)=0; 

Vzj=zeros(l,nt+l); 

Vyj=zeros(1,nt+1); 


%mass angular acceleration 
%mass angular velocity 
%mass angular displacement 
%pseudo-restoring force on mass 
%absolute horizontal position of mass 
%absolute vertical position of mass 
%global horizontal position of mass 
%global vertical position of mass 
%global horizontal position of support midpoint 
%global vertical position of support midpoint 
%pseudo-acceleration corresponding to the case 
break stability 

%global horizontal acceleration of support 

%global vertical acceleration of support 

%global horizontal velocity of support midpoint 
%global vertical velocity of support midpoint 


%linear horizontal velocity of mass during jump 
%linear vertical velocity of mass during jump 


%applying the data input in the support midpoint motion 


if dat==l 
Zs=data(: 
Ys=data(: 
else 

Zs=data(: 
Ys=zeros ( 
end 


,2) ' . /lOOO 
,3) ' . /lOOO 

,2) ' .71000 
1,nt+1); 
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%initializing the tension, support points position and string length 
%vectors 


Tright=zeros(1,nt+1); 
Tleft=zeros(l,nt+l); 
Zright=zeros(1,nt+1); 
Zleft=zeros(l,nt+l); 
Yright=zeros(1,nt+1); 
Yleft=zeros(l,nt+l); 
Lright=zeros(1,nt+1); 
Lleft=zeros(l,nt+l); 


%right string tension 
%left string tension 

%horizontal position of right support point 
%horizontal position of left support point 
%vertical position of right support point 
%vertical position of left support point 
%distance of mass from right support point 
%distance of mass from left support point 


for it=l:nt-l; %calculating the support displacements and 

%accelerations based on the support motion 
Zsdot (1, it + 1) = (Zs (1, it + 1) -Zs (1, it) ) /dt; 

Ysdot (1, it + 1) = (Ys (1, it + 1) -Ys (1, it) ) /dt; 

Zsddot(1,it + 1) = (Zsdot(1, it + 1)-Zsdot(1, it))/dt; 

Ysddot(1,it + 1) = (Ysdot(1,it + 1)-Ysdot(1, it))/dt; 

end 

vel(l,l)=0; %initial angulat velocity condition 

disp(l,l)=0; %initial angular displacement condition 

Zabs (1, 1) =-l* (sin (abs (disp (1,1) ) +a) -sin (a) ) ; 

Yabs(1, 1) =11-1*(cos(abs(disp(1,1))+a)); 

if abs(max(Ysddot))==0 && abs(max(Zsddot))==0&& abs(disp(1,1)==0) 
%if there is no acceleration applied to the support points 
%there is no delay for breaking the stability 
T=l; 

else 


for it=l:nt; 

acc_rest(1, it) = (abs(Zsddot(1,it)))/tan (a) -Ysddot(1,it); 

Tright(l,it) = (m/2)*((Ysddot(1,it)+g)/cos (a) +Zsddot(l,it)/sin(a)); 
Tleft(1,it)=(m/2)*((Ysddot(1,it)+g)/cos(a) -Zsddot(1,it)/sin(a)); 

Zright(l,it)=Zs(l,it)+l*tan(a) ; 

Zleft(1, it)=Zs(1,it)-l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1,it)=Ys (1,it); 

if acc_rest(1,it)<g %if the combined applied 

%acceleration is less than the required to break the stability 

disp(1,it)=0; %there is no oscillation 

vel (1,it)=0; 
acc(1,it)=0; 

Zabs (1, it) =-l* (sin (disp (1, 1) +a) -sin (a) ) ; 

Yabs(1, it) =11-1*(cos(abs(disp(1,1))+a)); 
Z(l,it)=Zabs(l,it)+Zs(l,it); 

Y(1,it)=Yabs(1,it)+Ys(1,it)-11; 

bright(1,it)=sqrt((Z(1,it)-Zright(l,it)).^2+(Y(l,it)- 
Yright (1, it) ) . ''2) ; 
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Lleft(1,it)=sqrt((Z(l,it)-Zleft(l,it)).^2+(Y(l,it)- 
Yleft (1, it) ) . ■'2) ; 

else 

T=it; 

Zabs (1, it) =-l* (sin (disp (1, 1) +a) -sin (a) ) ; 

Yabs(1, it) =11-1*(cos(abs(disp(1,1))+a)); 
Z(l,it)=Zabs(l,it)+Zs(l,it); 

Y (1, it) =Yabs (1, it) +Ys (1, it) -11; 

bright(1,it)=sqrt((Z(1,it)-Zright(l,it)).^2+(Y(l,it)- 
Yright (1, it) ) . ''2) ; 

Lleft (1, it) =sqrt ( (Z(l,it)-Zleft (l,it) ) .''2+(Y(l,it)- 
Yleft(1,it)) . "2); 

Vzl=Zsdot(1,T); 

Vyl=Ysdot (1,T); 

ACC=Zsddot(1,T); 

break 

end 


end 


end 


%Step after the stability is broken 


for it=T:nt; %loopl - initiating calculations after 

%stability is brocken 

if (bright(1, it)>0 | | Tleft(1,it)>0) %conditionl - oscillations take 
%place only if at least one of the tensions has a non-zero value 
force(1,it)=- 

(omega2+Ysddot (1, it) /I) * (sin (a) ) *cos (disp (1, it) ) *sign (disp (1, it) ) ; 
acc(1,it)=force (1,it)- 

(omega2+Ysddot(1, it)/I)*sin(disp(1, it))*cos(a) + (Zsddot(1,it))*((cos(abs 
(disp (1, it) ) +a) ) ) /I; 

vel(1, it + 1)=vel(1,it)+acc(1,it)*dt; 

disp(1, it + 1)=disp(1, it)+vel(1,it + 1)*dt; 

if disp(1,it+1)*disp(1,it)>=0 %condition2 - calculations 

%while the mass oscillates without changing strings 
disp(1, it + 1)=disp(1,it)+vel(1,it + 1)*dt; 

Zabs(l,it+1)=-sign((disp(1,it+1)))*1*(sin(abs(disp(l,it+l))+a)- 
sin (a)); 

Yabs(1,it+1)=11-1*(cos(abs(disp(1,it+1))+a)); 

Z(1, it + 1)=Zabs(1,it + 1)+Zs(1,it + 1); 

Y(1, it + 1)=Yabs(1,it + 1)+Ys(1,it + 1)-11; 

Zright(l,it)=Zs(l,it)+l*tan(a) ; 

Zleft(1, it)=Zs(1,it)-l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1, it)=Ys(1,it); 
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Lright(1,it)=sqrt((Z(1,it)-Zright(l,it)).^2+(Y(l,it)- 
Yright (1, it) ) . ''2) ; 

Lief t (1, it) =sqrt ( (Z(l,it)-Zleft (l,it) ) .''2+(Y(l,it)- 
Yleft(1,it)) . "2); 


else %condition2 - the mass reached 

%a string switch point 
Tl=it; 

V=vel(1,Tl); 

Vaft=sign(V)*sqrt((V*sin(pi/2 -2*a)) .^2+(r*V*cos (pi/2 -2*a)).^2); 
vel(1,Tl+1)=Vaft; 

disp(1, Tl + 1)=disp(1, it)+vel(1,Tl + 1)*dt; 

Zabs(l,it+l)=-sign((disp(l,it+l)))*1*(sin(abs(disp(l,it+l))+a)- 
sin (a)); 

Yabs(1,it+1)=11-1*(cos(abs(disp(1,it+1))+a)); 

Z(1, it + 1)=Zabs(1,it + 1)+Zs(1,it + 1); 

Y(1, it + 1)=Yabs(1,it + 1)+Ys(1,it + 1)-11; 

Zright(l,it)=Zs(l,it)+l*tan(a) ; 

Zleft(1, it)=Zs(1,it)-l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1, it)=Ys(1,it); 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (l,it)) .''2+(Y(l,it)- 
Yright(1,it)) . ^2); 

Lief t (1, it) =sqrt ( (Z(l,it)-Zleft (l,it) ) .''2+(Y(l,it)- 
Yleft(1,it)) . "2); 

end %condition2 

if disp(1,it+1)>0 %condition3 - calculating 

%the strings tension depending on which one of the strings is taught 

Tright(1, it + 1)=m*(vel(1, it + 1)*vel(1,it + 1)*l + Zsddot(1,it + 1)*sin(abs(disp 
(1, it + 1) ) +a) + (g+Ysddot(1,it + 1))*cos(abs (disp (1,it + 1)) +a) ) ; 

Tleft(1,it+1)=0; 
if Tright(1,it+1)>0 

Tright(1, it + 1)=Tright(1,it + 1); 
else 

Tright(1,it+1)=0; 

end 

else %condition3 

Tleft(1, it + 1)=m*(vel(1,it + 1)*vel(1,it + 1)*1- 
Zsddot(1,it+1)*sin(abs(disp(1,it+1))+a)+(g+Ysddot(1,it+1))*cos(abs(disp 
(l,it+l))+a)); 

Tright(1,it + 1)=0 ; 
if Tleft(1,it+1)>0 

Tleft(1,it + 1)=Tleft(l,it + l) ; 
else 

Tleft(1,it+1)=0; 

end 

end %condition3 

else %conditionl - the first 

%jump starts 
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T2=it 


V=vel(1,T2); 

disp(1, T2 + 1)=disp(1,T2)+vel(1,T2)*dt; 

Zabs(l,it)=-sign((disp(1,it)))*1*(sin(abs(disp(l,it))+a)-sin(a)); 
Yabs(1, it + 1)=11-1*(cos(abs(disp(1,it))+a)); 
Z(l,it)=Zabs(l,it)+Zs(l,it); 

Y (1, it) =Yabs (1, it) +Ys (1, it) -11; 

Zright(l,it)=Zs(l,it)+l*tan(a); 

Zleft (1, it) =Zs (1, it) -l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1,it)=Ys (1,it); 

Lright(1, it)=sqrt((Z(1,it)-Zright (l,it)) .^2+(Y(l,it)- 
Yright (1, it) ) . ''2) ; 

Lleft(1,it)=sqrt((Z(l,it)-Zleft(l,it)).^2+(Y(l,it)- 
Yleft (1, it) ) . ■'2) ; 

if disp(l,T2)>0 %condition4 - calculating 

%the linear velocity of the mass at the moment where the jump initiates 
Vmagn=abs(l*vel(1,T2)); 

Vz2=-sign(V)*Vmagn*cos((disp(1,T2)+a))+Zsdot(1,T2-1); 

Vy2=sign(V)*Vmagn*sin((disp(1,T2)+a))+Ysdot(1,T2-1); 
break 

else %condition4 

Vmagn=l*abs(vel(1, T2)) ; 

Vz2=-sign(V)*Vmagn*cos((abs(disp(1,T2))+a))+Zsdot(1,T2-1); 
Vy2=-sign(V)*Vmagn*sin((abs(disp(1,T2))+a))+Ysdot(1,T2-1); 
break 

end %condition4 

end %conditionl 


end %loopl 

%correcting the computation rounding errors 
LLright=Lright(l,T2)-0.004; 

LLleft=Lleft(1,T2)-0.004 ; 

Lright(1,T2)=Lright(l,T2)-0.004; 

Lleft(1, T2)=Lleft(1,T2)-0.004; 


while it<nt 
jump 

afterjump 
end 


time =data(:,1); 
if dat==l 
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figure (4) 

plot(time, disp,time,data(:,4), '— r') 
xlabel( 't (sec)' );ylabel( '\theta(rad)' ); 
axis([0 tf -pi/2 pi/2]) 
grid on 
figure (5) 

subplot(2, 1, 1);plot(time,Zabs,time,data(:,10) ./lOOO, '—r' ) 
xlabel ('t (sec) ');ylabel (' absolute horizontal position (m)' ); 
grid on 

axis ( [ 0 tf -1 1]) 

subplot(2,1,2);plot(time,Yabs,time,data(:,9)./1000+11 ,' —r' ) 
xlabel ('t (sec) ');ylabel (' absolute vertical position (m)' ); 
grid on 

axis ( [0 tf -1 1]) 
figure (6) 

subplot(2,1,1);plot(time,Z,time,data(:,8)./lOOO ,' —r' ) 
xlabel ('t (sec) '); ylabel (' total horizontal position (m) ' ); 
grid on 

axis ( [0 tf -1 3]) 

subplot(2,1,2);plot(time,Y,time, data(:, 7)./lOOO ,' —r' ) 
xlabel ('t (sec) ');ylabel (' total vertical position (m)' ); 
grid on 

axis([0 tf -1 1]) 
else 

figure (4) 

plot(time, disp,time,data(:,3), '— r') 
xlabel( 't (sec)' );ylabel( '\theta(rad)' ); 
axis([0 tf -pi/2 pi/2]) 
grid on 
figure (5) 

subplot(2,1,1);plot(time,Zabs,time,data(:,4) ./1000-Zs', '—r' ) 
xlabel ('t (sec) ');ylabel (' absolute horizontal position (m) ' ) ; 
grid on 

axis ( [0 tf -1 1]) 

subplot(2,1,2);plot(time,Yabs,time,data(:,5)./1000+ll-Ys' ,' —r' ) 
xlabel ('t (sec) '); ylabel (' absolute vertical position (m) ' ); 
grid on 

axis ( [0 tf -1 1]) 
figure (6) 

subplot(2,1,1);plot(time,Z,time,data(:,4)./1000 ,' —r' ) 
xlabel ('t (sec) ');ylabel (' total horizontal position (m)' ); 
grid on 

axis([0 tf -1 5] ) 

subplot(2,1,2);plot(time,Y,time, data(:, 5)./lOOO ,' —r' ) 
xlabel ('t (sec) ');ylabel (' total vertical position (m)' ); 
grid on 

axis ( [ 0 tf -1 0]) 
end 

figure (7) 

subplot(2,1,1);plot(time, bright ) 

xlabel ('t (sec) ');ylabel (' absolute horizontal position (m)' ); 
grid on 


126 


axis([0 tf 0 2 ] ) 

subplot(2,1,2);plot(time, Lleft ) 

xlabel('t (sec) '); ylabel (' absolute vertical position (m)' ); 
grid on 

axis ( [ 0 tf 0 2]) 


figure (8) 

subplot(2, 1 , 1 );plot(time,Zs ) 

xlabel('t (sec) '); ylabel( ' horizontal support position (m)' ); 
grid on 

axis ( [0 tf 0 5]) 

subplot(2,1,2);plot(time,Ys ) 

xlabel('t (sec) '); ylabel( 'vertical support position (m)' ); 
grid on 

axis ( [ 0 tf -1 2 ] ) 


figure (10) 

plot(time,Tright, time, Tleft) 

xlabel( 't (sec)' ); 

ylabel (' tension (N) ') 

legend (' right stringleft string') 

axis tight 


%jump.m subroutine - part of the Program_3.m code 


%Computing the response of the driven 2-pendulum when a jump is 
initiated 


"6 

"6 

%collecting data from the main file 

"6 

"6 

Lright(1,T2)=LLright; 

Lleft(1,T2)=LLleft; 

Vyj(l,T2)=Vy2; 

Vzj (1, T2)=Vz2; 
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%loopl - initiating 


for it=T2+l:nt; 

%calculations 

if (Lieft(1,it-1)>1 I I Lright(1,it-1)>1) %conditionl - if the 

%length of at ieast one of the strings is equal to 1, the the jump is 
%terminated 

T3=it; 

Zright(l,it)=Zs(l,it)+l*tan(a) ; 

Zleft (1, it) =Zs (1, it) -l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1, it)=Ys(1,it); 

Vz j(1,it)=Vz2; 

Vyj(1, it)=Vyj(1,it-1)-(g)*dt; 

Z(1, it) =Z(1, it-1)+Vz2*dt; 

Y (1, it) =Y (1, it-1) +Vyj (1, it) *dt; 

Zabs(l,it)=Z(l,it)-Zs(l,it); 

Yabs(l,it)=ll+Y(l,it)-Ys(l,it); 

vel(1,it)=0; 

disp(1,it)=0; 

Vz jl=Vz j(1,it); 

Vyjl=Vyj(1,it); 

Lright(1, it)=sqrt((Z(1,it)-Zright (l,it)) .^2+(Y(l,it)- 
Yright (1, it) ) . ''2) -0.004; 

Lleft(1,it)=sqrt((Z(l,it)-Zleft(l,it)).^2+(Y(l,it)- 
Yleft (1, it) ) . ■'2) -0.004; 

Vmagnj=sqrt(Vzjl*Vzjl+Vyjl+Vyjl) ; 

break 

else 

Vzj(1,it)=Vz2; %conditionl - the mass 

%follows a projectile path 

Vyj(1,it)=Vyj(1,it-1)-(g)*dt; %during the jump 

Z(1, it) =Z(1,it-1)+Vz2*dt; 

Y(l,it)=Y(l,it-l)tVyj(1,it)*dt; 

Zabs(l,it)=Z(l,it)-Zs(l,it); 

Yabs(l,it)=ll+Y(l,it)-Ys(l,it); 

Zright(l,it)=Zs(l,it)+l*tan(a) ; 

Zleft (1, it) =Zs (1, it) -l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1,it)=Ys (1,it); 

Lright(1, it)=sqrt((Z(1,it)-Zright (l,it)) .^2+(Y(l,it)- 
Yright (1, it) ) . ''2) -0.004; 

Lleft (1, it) =sqrt ( (Z(l,it)-Zleft (l,it) ) .''2+(Y(l,it)- 
Yleft (1,it)) ."2)-0.004; 

vel(1,it)=0; 

disp(1,it)=0; 

end %conditionl 

end %loopl 


if Lright(1,T3)>Lleft(1,T3) 

%condition2 - calculating the angular displacement 
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phi=atan((Zright(1,T3)-Z(1,T3))/(Yright(1,T3)-Y(1,T3))); %and 

%the velocity of the mass after the 

alpha=atan(Vyjl/Vzjl); %jump 

%is terminated depending on which string is taught 

vel(1,T3)=(1/1)*(-Vzjl*cos(phi)+Vyjl*sin(phi)); 
disp(1,T3)=phi-a; 

Tright(1,T3)=m*(vel(1,T3)*vel(1,T3)*l+Zsddot(l,T3)*sin(abs(disp(l,T3))+ 
a) + (g+Ysddot(1,T3))*cos(abs(disp(l,T3))+a)) ; 

Tleft(1,T3)=0; 

if Tright(1,it+1)>0 %condition3 - ensuring 

%that the string doesn't support compressive loads 
Tright(1,it+1)=Tright(1,it+1); 

else 

Tright(1,it + 1)=0 ; 

end 

else %condition2 

phi=atan((Z(1,T3)-Zleft(1,T3))/(Yleft(1,T3)-Y(1,T3))); 
phit=-phi; 

alpha=atan(Vzjl/Vyjl) ; 
alpha*180/pi; 

vel (1,T3) = (1/1)*(-Vzjl*cos(phi)-Vyjl*sin(phi) ) ; 
disp (1,T3)=-(phi-a); 

Tleft(1,T3)=m*(vel(1,T3)*vel(1,T3)*l+Zsddot(l,T3)*sin(abs(disp(l,T3))+a 
) + (g+Ysddot(1,T3))*cos(abs(disp(l,T3))+a)) ; 

Tright(1,T3)=0; 

if Tleft(1,it+1)>0 %condition3 

Tleft(1,it+1)=Tleft(1,it+1); 

else 

Tleft(1,it+1)=0; 

end 

end %condition2 
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%afterjump.m subroutine - part of the Program_3.m code 


%Computing the response of the driven 2-pendulum after a jump is 
%terminated 


for it=T3:nt; %loopl - initiating calculations after stability is 
%broken 

if (Tright(1, it)>0 | | Tleft(1,it)>0) %conditionl - oscillations take 
%place only if at least one of the tensions has a non-zero value 

force(1,it)=- 

(omega2+Ysddot (1, it) /I) * (sin (a) ) *cos (disp (1, it) ) *sign (disp (1, it) ) ; 
acc(1, it)=force (1,it)- 

(omega2+Ysddot (1, it) /I) *sin (disp (1, it) ) *cos (a) + (Zsddot (1, it) ) * ( (cos (abs 
(disp (1, it) ) +a) ) ) /I; 

vel(1, it + 1)=vel(1,it)+acc(1,it)*dt; 

disp(1, it + 1)=disp(1, it)+vel(1,it + 1)*dt; 

if disp(1,it+1)*disp(1,it)>=0 %condition2 - calculations 

%while the mass oscillates without changing strings 

disp(1, it + 1)=disp(1,it)+vel(1,it + 1)*dt; 

Zabs(l,it+1)=-sign((disp(1,it+1)))*1*(sin(abs(disp(l,it+l))+a)- 
sin (a)); 

Yabs(1,it+1)=11-1*(cos(abs(disp(1,it+1))+a)); 

Z(1, it + 1)=Zabs(1,it + 1)+Zs(1,it + 1); 

Y(1, it + 1)=Yabs(1,it + 1)+Ys(1,it + 1)-11; 

Zright(l,it)=Zs(l,it)+l*tan(a) ; 

Zleft (1, it) =Zs (1, it) -l*tan (a) ; 

Yright(1,it)=Ys(1, it) ; 

Yleft(1, it)=Ys(1,it); 

Lright (1, it) =sqrt ( (Z (1, it) -Zright (l,it)) .''2+(Y(l,it)- 
Yright(1,it)) . ^2); 

Lleft(1, it)=sqrt((Z(l,it)-Zleft (l,it)) .^2+(Y(l,it)- 
Yleft (1, it) ) . ■'2) ; 


else %condition2 - the mass reached 

%a string switch point 

T2=it; 

V=vel(1,T2); 

Vaft=sign(V)*sqrt((V*sin(pi/2 -2*a)) .^2+(r*V*cos (pi/2 -2*a)).^2); 
vel(1,T2+1)=Vaft; 

disp(1, T2 + 1)=disp(1, T2)+vel(1,T2 + 1)*dt; 

Zabs(l,it+1)=-sign((disp(1,it+1)))*1*(sin(abs(disp(l,it+l))+a)- 
sin (a)); 
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Yabs(1,it+1)=11-1*(cos(abs(disp(1,it+1))+a)); 

Z (1,it + 1)=Zabs(1,it + 1)+Zs(1,it + 1); 

Y(1,it+1)=Yabs(1,it+1)+Ys(1,it+1)-11; 

Zright(l,it)=Zs(l,it)+l*tan(a); 

Zleft(1, it)=Zs(1,it)-l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1, it)=Ys(1,it); 

bright (1, it) =sqrt ( (Z (1, it) -Zright (l,it)) .''2+(Y(l,it)- 
Yright(1,it)) . ^2); 

Lief t (1, it) =sqrt ( (Z(l,it)-Zleft (l,it) ) .''2+(Y(l,it)- 
Yleft (1,it)) . "2) ; 


end %condition2 

if disp(1,it+1)>0 %condition3 - calculating 

the strings tension depending on which one of the strings is taught 

Tright(1, it + 1)=m*(vel(1, it + 1)*vel(1,it + 1)*l + Zsddot(1,it + 1)*sin(abs(disp 
(1, it + 1) ) +a) + (g+Ysddot (1, it + 1) ) *cos (abs (disp (1, it + 1) ) +a) ) ; 

Tleft(1,it+1)=0; 

if Tright(1,it+1)>0 %condition4 - ensuring that 

%the string doesn't support compressive loads 

Tright(1,it+1)=Tright(1,it+1); 
else 

Tright(1,it + 1)=0 ; 

end 

else %condition3 

Tleft(1, it + 1)=m*(vel(1,it + 1)*vel(1,it + 1)*1- 
Zsddot(1,it+1)*sin(abs(disp(1,it+1))+a)+(g+Ysddot(1,it+1))*cos(abs(disp 
(l,it+l))+a)); 

Tright(1,it + 1)=0 ; 

if Tleft(1,it+1)>0 %condition4 - ensuring that 

%the string doesn't support compressive loads 

Tleft(1,it + 1)=Tleft(l,it + l) ; 
else 

Tleft(1,it+1)=0; 

end 

end %condition3 

else %conditionl 

T2=it 

V=vel(1,T2); 

disp(1, T2 + 1)=disp(1,it)+vel(1,T2)*dt; 

Zabs(l,it)=-sign((disp(1,it)))*1*(sin(abs(disp(l,it))+a)-sin(a)); 
Yabs(1, it + 1)=11-1*(cos(abs(disp(1,it))+a)); 
Z(l,it)=Zabs(l,it)+Zs(l,it); 
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Y (1, it) =Yabs (1, it) +Ys (1, it) -11; 

Zright(l,it)=Zs(l,it)+l*tan(a); 

Zleft(1, it)=Zs(1,it)-l*tan (a) ; 

Yright(l,it)=Ys(l,it); 

Yleft(1,it)=Ys (1,it); 

bright(1, it)=sqrt((Z(1,it)-Zright (l,it)) .^2+(Y(l,it)- 
Yright (1, it) ) . ''2) ; 

Lleft (1, it) =sqrt ( (Z(l,it)-Zleft (l,it) ) .''2+(Y(l,it)- 
Yleft(1,it)) . "2); 

if disp(1,T2+1)>0 %condition5 - calculating 

%the velocity when ont of the strings becomes taught 
Vmagn=abs(l*vel(1,T2)); 

Vz2=-sign(V)*Vmagn*cos((disp(1,T2)+a))+Zsdot(1,T2-1) ; 

Vy2=sign(V)*Vmagn*sin((disp(1,T2)+a))+Ysdot(1,T2-1); 
break 

else %condition5 

Vmagn=abs(l*vel(1, T2 ) ) ; 

Vz2=-sign(V)*Vmagn*cos((abs(disp(1,T2))+a))+Zsdot(1,T2-1); 
Vy2=-sign(V)*Vmagn*sin((abs(disp(1,T2))+a))+Ysdot(1,T2-1); 
break 


end 

end 


%condition5 

%conditionl 


end 


%loopl 


vel(1,T2); 

Vz2; 

Vy2; 

LLright=Lright(l,T2)-0.004; 
LLleft=Lleft(1, T2)-0.004 ; 
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